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Chapter 1 
INTRODUCTION 

This report augments the annual report for 1981 that des- 
cribes the work performed by the Virginia Tech Satellite 
Communications Group under Jet Propulsion Laboratory Con- 
tract No. 955954. This report details the modeling phase of 
the overall effort. A separate report was necessary because 
of the already excessive length of the annual report. 

This year's effort in modeling consisted of five tasks. 
These tasks and a brief report on the progress made in each 
follow. 

1 . Simple Attenuation Model . The development of the 
simple attenuation model is complete. It consists 
of an exponentially shaped spatial rain rate dis- 
tribution. Many comparisons to measured data and 
to other models have been made. See Chapter 3 for 
a more complete discussion. 

2. Attenuation Exceedance Testi ng. This effort is 
the prediction of rain attenuation exceedance sta- 
tistics. The procedure is to couple a rain rate 
exceedance model (such as that recommended by 
CCIR) with the simple attenuation model which 


predicts slant path attenuation for a given point 
rain rate. Results of this technique are also 
found in Chapter 3. 


Isolation versus Attenuation . After the spatial 
rain rate distribution is established it can be 
used together with a complete depolarization com- 
putational model to calculate isolation (as a 
function of attenuation). This particular task is 
incomplete because of the extensive effort 
required on the development of the spatial rain 
rate distribution (Task A) and the development of 
the multiple scattering model (Task E), both of 
which must precede this task. 

Single Particle Scattering Computations . In order 
to make complete depolarization calculations for 
rain media the single rain drop scattering coeffi- 
cients must be evaluated at the frequency of 
interest. This task was completed by developing 
and testing a computer program for calculation of 
oblate spheroidal rain drops. It uses the Fre- 
dholm integral equation method. Unfortunately, 
the calculations are very complicated and involve 
considerable computer time. Complete details are 
presented in Chapter 4. 
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5. 


Multiple Scattering Model . The derivation and 
initial testing of a general rain depolarize mi on 
model which includes multiple scattering effects 
has been completed under this task. Chapter 5 
discusses this effort in detail. In Chapter 2 the 
multiple scattering model is placed into perspec- 
tive relative to other computing capabilities. 



Chapter 2 

MODELING CAPABILITIES 

In the past several years three levels of rain propaga- 
tion computing capabilities have been developed. It may be 
helpful to summarize this capability before proceeding into 
the details of the recent findings. 

1. SAM - Simple Attenuation Model. This model is 
intended for use in predicting rain attenuation 
encountered on earth-space communication links. 

Given the earth terminal location, attenuation 
versus point rain rate is easily calculated. 

Further, input of the rain rate exceedance (such 
as the appropriate global rain rate region adopted 
by CCIR) permits prediction of attenuation excee- 
dance (percent time during the year a given atten- 
uation level is exceeded) on a worldwide basis. 
Although computing capability over the frequency 
range of 1 to 1000 GHz exists, only the 11 to 35 
GHz band has been tested. 

2. KPP - Rain Propagation Prediction Model. This is 
a first-order multiple scattering model developed 
for computation of isolation and phase as well as 
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attenuation as a function of point rain fall rate 
on an earth-space link. The program can be used 
for instantaneous or average computations by 
entering the appropriate spatial rain rate profile 
in a piecewise (ten level) uniform manner. The 
restriction on the capacity of the method is that 
computations are restricted to frequencies for 
which single rain drop scattering coefficients are 
available. Currently the program operates at fre- 
quencies of 11, 14, 20, and 30 GHz. 

3. Multiple Scattering Model. This program operates 
in essentially the same manner as RPP with the 
same inputs and outputs. The algorithm is, how- 
ever, more general in that all orders of multiple 
scattering are included. The recent computer pro- 
gram version of this technique has shown that num- 
ercial results are very similar to those of the 
RPP first-order multiple scattering program for 
frequencies of 30 GHz or under. Differences are 
expected to occur for higher frequencies. The 
multiple scattering computer program is not only 
more general than RPP but is also computationally 
more efficient. 

The modeling capabilities are illustrated in the block 
diagram of Pig. 1. Since the multiple scattering model is 
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Input: 

SAM 

Link characteristics and location 

Output : 

Rain rate exceedence 
(for %1 vs A Prediction) 
A vs R, T % vs A 

Restrictions: 1 to 1000 GHz 


(Tested in 11-35 GHz range) 


First-order multiple scattering with cellular input 

Input: Link characteristics 

Rain profile-piecewise uniform 

Output: Rj Ai I, $ 

Restrictions: 11, 14, 20, 30 GHz 


Multiple Scattering Program 


Operates in essentially same manner as 
RPP, but is more efficient and includes 
multiple scattering effects. 


Single Drop Scattering Coefficient Program | 


Run for several drop sizes at a given 
frequency. Curve fit to give coefficient 
as function of drop size. 


Figure 1. Summary of rain propagation modeling capabilities at Virginia 













Chapter 3 

THE SIMPLE ATTENUATION MODEL 


This chapter is similar to the paper in the Nov/Dec 1982 
special issue of Radio Science on the NASA Propagation Pro- 
gram. More details can be found in the interim report 
"Estimation of Rain Attenuation on Earth-Space Millimeter 
Wave Communication Links," by W. K. Dishnun and W. L. Stutz- 
man published May 1982. 


3.1 INTRODUCTION 

As earth-satellite communications increase, economic con- 
siderations become more important. One method of reducing 
system cost is to operate with lower signal power margin. 
Accurate calculation of predicted signal power budgets per- 
mit systems to operate with a narrower margin for fading. 
Central tc such calculations for links operating above 10 
GHz is the accurate modeling of rain fading [Crane, 1977; 
Ippolito, 1981]. Initial attenuation prediction attempts 
involved extrapolation of measurements to other locations, 
frequencies, and elevation angles. The complex nature and 
regional variability of rain make this approach highly inac- 
curate. Over the past several years research activity has 
been very vigorous with many models being proposed in an 


attempt to improve predictions, A wealth of literature 
exists and the reader is referred to several review papers 
[Rogers, 1976; Crane, 1977; Lane and Stutsman, 1980b; Brus- 
saard, 1981; Ippolito et al . , 1982], 

There are some areas of attenuation modeling which are 
incomplete due to the lack of sufficient physical data. 
There is also some disagreement among researchers about how 
the problem should be attacked. These not withstanding, 
research investigations over the past several years have, 
indeed, moved closer together in approach. As pointed out 
by Brussard [1981], rain propagation research has the goal 
of providing information useful for communication system 
design. With this is mind the following have been identi- 
fied by Fedi [1981a] as being desirable features of a pred- 
iction method: 

1 . Simple . The model should be easy to apply to com- 
munication system calculations. The unnecessary 
introduction of new parameters and mathematical 
complexity is to be avoided. 

2. Physically sound . As much as possible, the model 
should be checked against directly observed physi- 
cal data, such as spatial rain behavior. 

3. Data Tested . The model should be tested against 
measured data from many different regions. Empha- 
sis should be given to the data at low percentages 


of time that are of most interest to system 
designers . 

4. Flexible . As more data becomes available and a 
deeper understanding is obtained, model refine- 
ments surely follow. The model should be struc- 
tured to accept modifications. 

Herein a prediction method that incorporates all of the 
features mentioned above is presented. In spite of converg- 
ing thought, it would be presumptuous to report that this 
work represents the prevailing trends of all researchers. 
However, even though differences are present in current 
models, the opportunity is taken to include many commonly 
accepted elements into a single model. Complete details on 
this investigation are given by Dishman and Stut 2 man [1982] 

In Section 3.2 of this report fundamental concepts of 
attenuation modeling are discussed. In Section 3.3 we deal 
with the difficult problem of describing the spatial distri- 
bution of rain and we propose an exponential rain rate pro- 
file. The complete attenuation model is presented in Sec- 
tion 3.4 and it is evaluated by comparison to measured data 
from many experiments around the world and to predictions 


from other models 



3.2 FUNDAMENTAL CONCEPTS IN ATTENUATION MODELING 


All attenuation models are "semi-empirical" in nature in 
that they employ attenuation data in the development of the 
model. This is done, however, to various degrees and we 
shall classify attenuation models as using either an empiri- 
cal approach, a rain-cell approach, or a rain profile 
approach. The empirical approach develops an expression for 
attenuation directly from measured attenuation. Models 
based on a completely empirical approach are usually easy to 
apply, but do not relate directly to rain parameters. The 
methods of Lin [1979] and the CCIR [1981b] are examples of 
empirical models. Models based on a rain cell approach 
include the physically realistic idea of a randomly located 
"cell". The rain-cell models of Misme and Waldteufel [1980] 
and Lane and Stutzman [1980a, 1980b] require computer pro- 
grams for evaluation. Models based on a rain profile employ 
an effective rain rate spatial distribution and are gener- 
ally easy to use. 

The decision of which approach to select for attenuation 
prediction is guided by the desirable features of a model 
discussed in Section 3.1. In particular, the model must be 
simple and physically sound. Empirical models are usually 
easy to use, but lack the necessary physical foundations, 
which in turn, leads to questions concerning the range of 
applicability. On the other hand, rain-cell models are more 
physically acceptable, but are generally more complex. It 
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has been found [Lane and Stutzman, 1980a, 1980b] that a 
rain-cell model which is stochastic in nature, allowing for 
the rain cell position to be random, is unnecessary and that 
a space-fixed rain profile is sufficient for prediction of 
average attenuation. Thus, we select the rain profile model 
approach, which contains both the elements of being physi- 
cally sound and simple. 

Now, the rain attenuation modeling problem using the rain 
profile approach can be divided into three different areas 
[Fedi, 1981a]: (1) The relationship between specific atten- 

uation and rain rate; (2) The statistics of point rainfall 
intensity; and (3) The spatial distribution of rainfall. 

The first two areas are relatively well understood and have 
received much attention in the literature. A brief discus- 
sion of these topics will be presented in this section along 
with a discussion of the path integral concept. The spatial 
distribution of rainfall will be treated separately in the 
next section. 


3.2.1 Specific Attenuation 



been given by Olsen et al.[1978]. The main parameters asso- 
dated with the microstructure of rain are the shape, size 
distribution, and temperature of the raindrops. 

Raindrops are usually assumed to be either spherical, 
oblate spheroidal, or of the shape described by Pruppacher 
and Pitter [1971]. Use of the latter two shapes allows one 
to include the effects of wave polarization, drop canting 
angle, and slant-path elevation angle in the calculation of 
specific attenuation. Attenuation values computed assuming 
spherical drops generally lie between the extremes of values 
computed for vertical and horizontal linear polarization 
assuming distorted drop shapes. Because the errors between 
attenuation values computed using spherical drops and those 
computed using distorted drop shapes are typically 10% or 
less in the frequency range of interest [Crane, 1977; Olsen 
et al. , 1978], specific attenuation values calculated for 
spherical drops will produce adequate results. 

Various drop size distributions have been considered in 
the. calculation of specific attenuation. These include the 
distributions of Laws and Parsons [1943], Marshal and Palmer 
[1948], and Joss et al.[1968]. While there is little dif- 
ference in the values of actenuation computed using these 
distributions for frequencies below 30 GHz, the Laws and 
Parsons distribution is preferred because of the tendency of 
the other distributions to overestimate attenuation at the 
higher frequencies [Crane, 1977? Olsen et al-, 1978; Upton 
et al., 1980; Ippolito, 1981]. 


The temperature of the raindrops is perhaps the most 
critical parameter. The assumed temperature has little 
effect above 15 GHz, but specific attenuation values are 
very sensitive to temperature variations in the 11-14 GHz 
band [Olsen et al . , 1978; Upton et al . , 1980? Thompson et 
al. , 1980]. Values of a and b are generally available for 
temperatures of 0° C and 20° C. While 20° C is a reasonable 
assumption for terrestrial link attenuation prediction 
[Damosso et al . , 19802, it is probably not representative of 
the temperature for earth-space paths [Thompson et al - , 

1980]. Attempts have been made to include the temperature 
variation of specific attenuation with altitude [Misme and 
Waldteuf ul , 1980], but this introduces a complexity into the 
calculation of attenuation that has little effect upon the 
results. For most climates, the assumption of 0° C drop 
temperature should give good results [Olsen et al., 1978] „ 
Based upon the previous assumptions about the fine struc- 
ture of rainfall, we believe that the use of spherical 
drops, the Laws and Parsons drop size distribution, and 0° C 
rain temperature will give reasonably accurate values of 
specific attenuation. A convenient source of computing the 
coefficients a and b for any frequency of interest are the 


.-FT* U .*H' 


ORIGINAL PAGE 19 
OF POOR QUALITY 


following equations taken from Olsen et al . [ 1978 ] : 


a(f) - 


4.21 x 10~ 5 f 2 * 42 


4.09 x io' 2 f 0 * 699 


2.9 < f < 54 GHz 


54 < f < 180 GHz 


( 2 ) 


b(f) = 


1.41 f 


2.63 f 


-0.0779 


-0.272 


8.5 < f < 25 GHz 


25 «£ f < 164 GHz 


C3) 


Should more precise values of a and b be desired, the 
tabulated values given in [Olsen et al., 1978] should be 
used. To include the effects of polarization and elevation 
angle, refer to the tabulated values recommended by the CCIR 
[1981c] or to the regression equations given by Thompson et 
al. [1980] and Damosso [1981]. 

3.2.2 Point Rainfall Intensity Distribution 

h critical parameter in the estimation of attenuation 
exceedance is the point rain rate distribution. As pointed 
ou' by Crane [1977], cumulative rain rate distributions may 
show considerable variability from year to year. For this 
reason care must be taken when estimating the average dis- 
tribution at a site. An excellent review of procedures for 
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PERCENT TIME RAIN RATE IS EXCEEDED 


imUtflAi 


I f I 


• • • • Measured Data 

Predicted Using 

COIR Region K 


RAIN RATE (mm/hr) 

Cumulative distribution of rainfall intensity at Blacksburg , VA 
as Treasured for the period June 1976 to June 1979 (dots) and as 
estimated for rain ciimate zone K using the CCIR distribution 
[CCIR, 1981a]. 
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estimating the rain intensity distribution is given by Fedi 
[1981a]. For use in attenuation estimates, data obtained 
from local sources are preferred. However, when adequate 
local data is not available, the distributions can be esti- 
mated from the rain climate region maps recommended by the 
CCIR [1981a], These maps present cumulative rainfall dis- 
tributions for 14 different regions of the world. A compar- 
ison between three years of measured data from Blacksburg, 
VA, and the corresponding CCIR region is shown in Fig. 1. 

The measurements were taken with a tipping bucket type rain 
gauge which usually provides a good approximation to the 
instantaneous rain rate. In general we have obtained good 
results using the CCIR rain region model [Dishman and Stutz- 
man, 1982]. 


3.2.3 The Path Integral 

If the rain rate profile, R(l), were known along the 
extent of the propagation path, L, it would be a simple mat- 
ter to calculate the total attenuation by integrating over 
the incremental (or specific) attenuation: 


A(R 0 ) 


(L 

a[R(£)] d£ 

J 0 


(4) 


In this equation we have indicated that total attenuation is 
a function of the point rain rate at the 9,-0 end of the 
path, or R 0 = R(&=0). This is done to emphasize the real- 
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m 


ities of the problem. Rarely is R(M known, whereas point 


rain rate is directly measured and is available in histori- 


cal data form. Thus, A is expressed in terms of R 0 . As 


mentioned earlier, one step in the modeling process is to 


obtain a spatial rain rate profile. Once this is done (4) 


can be used to make predictions. Since the profile R( i) is 


an average one, the calculated attenuation will be also, 


Several methods have been proposed to utilize the path 


integral, but without explicitly developing a rain profile 


R(ji). These methods include using an effective path length 


L e (where A = a (R 0 ) L e ) , a path averaged rain rate 
(where A = a R ave L), or a path average factor r = 
R ave /R 0 . Although good results have been obtained in 


some cases, each of these approaches involve some approxima- 


tion which limits its usefulness CKheirallah et al . , 1980? 


Dishman and Stutzman, 1982]. Attenuation prediction for 


arbitrary situations follows directly from the path integral 


(4) and this is discussed in the next section. 


3.3 SPATIAL RAINFALL DISTRIBUTION 


The most difficult parameter of an attenuation model to 


characterize is the spatial distribution of rain. Gener- 


ally, precipitation systems are combinations of both strati- 


form and convective rain structures. Radar measurements 


indicate that most precipitation is characterized by large 


areas of low rates with a number of smaller regions of high 
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rain rates [Crane, 1977], It is the presence of these 
imbedded rain “cells* 1 that makes the spatial distribution of 
rain difficult to describe. Because rain rate statistics 
are usually only available for point rain rates, it is 
necessary to develop an "effective" spatial rain distribu- 
tion model (or rain profile) that relates the rainfall along 
a path to the rainfall at a point. A rain profile is not 
applicable to single event analysis? it is useful for commu- 
nication system design situations that involve long term 
performance. In other words, the wide variations in rain 
cells (observed for short time periods) tend to average out 
over long time periods, and a rain profile is useful in sta- 
tistical predictions. 

The rain profile for use in attenuation prediction must 
include both the horizontal and vertical spatial variations 
of rain as discussed in this section. Ideally, spatial var- 
iations should be determined with direct measurements of 
rain behavior, such as with rain gauge networks and radar. 
But there are not enough direct observations available to 
completely develop a model. The large indirect-measurement 
data base of attenuation on earth-space links must also be 
used [Fedi , 1981b]. Direct measurements are used in this 
section to establish a rain cell model which is exponential 
in shape. In the next section it is shown that the results 
agree with inferences from indirect measurements as well. 
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We will adopt the customary assumption that this effec- 
tive distribution is the same in all geographic regions of 
interest. It is noted, however, that in regions where oro 


graphic features play a strong role or where rainfall is 
dominated by a particular form of precipitation (e.g. typ- 
hoons, hurricanes) this assumption may not hold. 

3.3.1 Vertical Variation of Rain 

Radar observations have shown that the vertical structure 
of precipitation is characterized by two different regions. 
The upper region consists of a mixture of ice and snow and 
does not contribute significantly to attenuation at frequen- 
cies below 60 GHz [CCIR, 1981a]. The lower region is mostly 
rain and is the primary source of attenuation. The transi- 
tion height between the two regions corresponds approxi- 
mately to the height of the 0® C isotherm. 

Goldhirsh and Katz [1979] have presented median reflec- 
tivity factor profiles obtained from radar observations of 
summer rain cells at Wallops Island, VA. These profiles 
indicate that reflectivity is essentially constant up to a 
certain height and drops off rapidly above this height. 
Similar observations have been reported by other researchers 
[CCIR, 1981a]. This leads to the assumption of uniform rain 
structure from the ground up to an u e£ f ective" rain height 
H e , as shown in Fig. 2. The contribution to attenuation 


i 



I 



by particles above the effective rain height will be 
neglected. 






The radar reflectivity profiles given by Goldhirsh and 
Katz [1979] show that the rain height is approximately con- 
stant and equal to the height of the 0° C isotherm for low 
rain rates. This is consistent with other observations of 
stratiform rain. As rain -rate increases, however, the rain 
height indicated by the reflectivity profiles also 
increases. This increase is due to the structure of convec- 
tive rain cells in which liquid water may be carried well 
above the 0° C isotherm level by updrafts. The data of Gol- 
dhirsh and Katz indicates that on the average the rain 
height may extend approximately 1 km above the 0° C isotherm 
height for rain rates in excess of lOOmm/hr. A reasonable 
model for the effective rain height therefore consists of 
using the 0° C isotherm height for low rain rates and adding 
a rain rate dependent term to the 0° C isotherm height for 
higher rain rates. We propose the following simple rela- 
tionship for the "effective" rain height H e in km: 

f 

H. R < 10 mm/hr 

i o 


= 


H i + log 


R 

o 

10 


(5) 


R q > 10 mm/hr 


where R Q is the point rain rate in mm/hr and is the 0° 
C isotherm height in km. The breakpoint of 10 mm/hr was 
chosen because it corresponds to the approximate value of 
the maximum rain rate associated with stratiform rain. 
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The 0° C isotherm height varies with latitude and 
with the season of the year. Zonally averaged values of 
versus latitude for the four seasons are given by Oort 
and Rasmussen [19711. Based upon this data, Crane [1978] 
approximated the average height of H-^ by 


4.8 


A | < 30° 


H i = 


( 6 ) 


7.8 - 0 . 1 j A 


A | > 30° 


where A is the latitude in degrees. As stated by Crane, 
this expression is an approximation to the observed mean 
seasonal values that lies midway between the summer and 
spring or fall values. The effective rain height H e is 
shown in Fig. 3 as a function of latitude and rain rate. 
These curves were obtained using (6) in (5). 

3.3.2 Horizontal Variation of Rain 

Convective cells imbedded in stratiform rain render the 
distribution of rain nonuniform in the horizontal direction. 
Direct methods of observing rain cell structure include rain 
gauge networks and radar, and possibly the synthetic storm 
method. Radar observations yield the best information about 
the rain structure, but have not yet been used widely in 
determining the path-averaged to point rainfall relation- 
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LATITUDE (degrees) 

Figure 3. Variation of the effective rain height H e with latitude and 
point rainfall rate, R q . 
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ship. There have been a limited number of rain gauge net- 
works operated around the world. Sims and Jones [1975] ana- 
lysed data obtained from the 1946 Thunderstorm Project with 
lines of rain gauges in Florida and from a network of gauges 
in Illinois during the summer of 1970. Freeny and Gabbe 
[1969] reported on the results of a dense rain gauge network 
operated for six months in New Jersey. Harden et al . [1977] 
and Valentin [1977] operated rain gauge networks in conjunc- 
tion with radio links in England and West Germany, respec- 
tively. The synthetic storm method, as employed by Drufuca 
[1974] and others [Harden et al., 1974; Kheirallah et al., 
1980; Watson et al. , 1977] uses the translational velocity 
of a storm over a rain gauge to convert point rainfall sta- 
tistics to spatial distributions. Although this method is 
somewhat difficult to apply at an arbitrary site, it does 
provide useful information for developing a general rela- 
tionship for the spatial rain distribution. 

The point-to-path rainfall relation can be indirectly 
characterised from terrestrial link attenuation data. The 
most widely used method is the effective path length method, 
which was mentioned in Section 3.2.3. The exclusive use of 
these methods for determining spatial characteristics of 
rain should be avoided because of the non-linear relation- 
ship between specific attenuation and rain rate. However, 
when used in conjunction with data from rain gauge networks, 
attenuation data is very valuable. 
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The data of Sims and Jones [19753 and of Harden et 
al.[1977] both suggest that the point and path-averaged rain 
rates are the same up to rates of 10 to 14 mm/hr. Similar 
results were reported by Kheirallah et al.[1980] based upon 
synthetic storm studies at three locations in Canada. For 
rain rates in excess of 10 mm/hr, the path-averaged rain 
rate decreases as point rain rates increase and as path 
lengths increase. Several researchers [Valentin, 1977? Bar- 
sis and Samson, 1976] have presented curves derived from 
rain gauge data that illustrate the point-to-path rain rela- 
tionship. 

As pointed out by Fedi [1981a] and Kheirallah et 
al.[1980], the use of path-averaged rain rate in the spe- 
cific attenuation relationship is not directly applicable 
due to the influence of the exponent in (1). Crane [1980] 
attempted to overcome this problem by fitting a power-law 
relation to the rain gauge data from Europe and the United 
States, then differentiating these data to obtain the path 
profile of the rain. The resulting data were fitted with a 
series of exponential functions that could be used in the 
path integral relation (4). Kheirallah et al.[1980] point 
out that the Crane model overestimates the path averaged 
rain rates for low values of point rain rate. 

The piecewise uniform path profile model [Persinger et 
al . , 1980] gives good results but has only two rain rate 
values allowed along the path. The piecewise exponential 
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path profile of the global model [Crane# 1980] also gives 
good results but is unnecessarily complicated and does not 
include uniform rain rates for low point rain rate values. 
Combining these ideas together with a goal of simplicity we 
propose the following exponential -shaped effective path pro 
file for rain rate: 


R R < 10 irnn/hr 

o o 


Hz) 


-y In (R /10) z 

R o e 


. ( 7 ) 


R q > 10 mm/hr 


&o is the point rainfall intensity, z is the horizontal 
distance along the path, and y is a parameter controlling 
the rate of decay of the profile. 

The path-averaged rain rate for a path of length D is 
found from (7) as 


R 


ave 


= i f D 
D Jn 


R(z)dz = 


R 


R. 


1 - e 


~y ln(R o /10)D 


o y In (R7/10)D 


o' 


R q < 10 mm/hr 


( 8 ) 


R > 10 mm/hr 
o 


This relation was compared to measured values of path rain 
rate obtained from the literature [Sims and Jones, 1975? and 


- 28 - 


ORIGINAL FAGS 

OF POOR QUALITY 

Valentin, 1977; Freeny and Gabbe, 1969; Harden et al . , 1977] 
and to values derived using the synthetic storm model 
[Kheirallah et al . , 1980] in order to establish the value of 
the parameter y . Th j i re : * ^ bounded by values of y bet- 
ween 1/10 and 1/3 v, r. 'S '2 giving a best fit. A plot of 
the normalized ] f of (7) versus distance is shown 


in Fig. 4 for y 


v.Yf/1 several values of rain rate. 


3.3.3 Summarv of the Proposed Rain Rate Profile 


Equations (5)- (7) describe the proposed spatial distribu- 
tions in the vertical and horizontal planes. Because the 
rain is assumed to be uniform in the vertical direction up 
to H 0 , the rain profile R(M along the slant path (see 
Fig. 2) can be derived using the simple trigometric rela- 
tionship z = £ cos e in (7), giving 


■0) = 


-Y ln(R _/10)f. cos e 


R o e 


R < 10 mra/hr 
o 


R > 10 mm/hr 
o 


for l ^ L, where 


H - H 
e o 

sin e 


This expression is likely to be valid for elevation angles 
above about 10°. For lower elevation angles a value of L as 
suggested by the CCIR [ 1981b] could be used, although very 






Figure 4. 


Normalized rain p 
rainfall rate fro 


! 



^teggtega 

!l v 



distance and point 


ORJGflMA 



low elevation angle satellite data and terrestrial data have 
not been compared to this model. H 0 is the altitude of 
the earth station location and H e is defined by (5) and 
(6). This expression will be used in the next section to 
derive the corresponding total attenuation model. 


3.4 THE SIMPLE ATTENUATION MODEL AND ITS PERFORMANCE 
The total attenuation due to a point rainfall rate R 0 
is easily computed using the effective rain profile from (9) 
in the path integral (4). Evaluating the integral gives 


A(R 0 ) = < 



-Y b ln(R o /10)L cos e 
Y b ln(U Q /10) cos e 


R q < 10 mm/hr 


( 11 ) 


R Q > 10 mm/hr 


where the path length L is given by (10). 

The simple attenuation model (SAM) given by (11) is a 
function of the point rainfall intensity only; it is not a 
function of the percentage of time that rain rate is 
exceeded. This decoupling of the attenuation model and the 
rain rate statistics allows evaluation of the model parame- 
ters independent of the errors associated with the rain rate 
statistics . 
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Many propagation experiments using satellite beacons have 
been conducted in North America, Europe, and Japan. An 
extensive collection of the rain rate, attenuation, and iso- 
lation statistics associated with these experiments has been 
assembled at VPI&SU to form a data base for use in develop- 
ing and testing propagation models. The attenuation versus 
rain rate data (based on equal probability levels) from this 
data base were compared to the predictions calculated from 
(11). The best average agreement was found when y = 1/22. 
This is the same value found when comparing the path-aver- 
aged rain rate to directly measured rain gauge data (see 
Section 3.3.2). Thus, both direct and indirect data agree 
with the model of (11) when y = l/22 is used. This adds 
confidence that the spatial rain rate variation of Section 
3.3 is an accurate model. A plot of typical A vs R 0 data 
is shown in Fig. 5. The data were obtained from three years 
of measurements in Blacksburg, VA using the CTS (11.7 GHz) 
and COMSTAR (19 and 28 GHz) beacons. Also shown in Fig. 5 
is the attenuation calculated using (11) with y = 1/22. 

Estimates of attenuation statistics are found by combin- 
ing the attenuation versus rain rate model with a rainfall 
distribution (see Section 3.2.2). An example of this proce- 
dure is given in Fig. 6. The predicted attenuation distri- 
butions were calculated using the CCIR rain rate distribu- 
tion for the Blacksburg region (see Fig 4) together with the 
SAM predictions of attenuation as a function of rain rate. 
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Figure 5. Attenuation vs* rain rate data from the following experiments 
in Blacksburg, VA: CTS (11.7 GHz, 33 degree elevation angle, 

circular polarization, June 1976 to June 1979), COMSTAR (19.4 
and 28.56 GHz, 45 degree elevation angle, vertical polarization, 
July 1977 to August 1980) . The solid lines are the estimates for 
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To do this we make the customary assumption that the proba- 
bility of the attenuation exceeding a certain values is the 
same as the probability of the point rainfall intensity 
exceeding the point rainfall rate used to predict the atten- 
uation. Also plotted in Fig. 4 are the measured attenuation 
distributions obtained from three years of observations in 
Blacksburg . 

The complete VPI&SU data base which consists of attenua- 
tion measurements from 47 experiments is presented in Table 
1 [Dishman and Stutzman, 1982]. The experiments represent 
17 different sites in the U.S., Europe, and Japan ranging in 
latitude from 28 to 52° N, varying in frequency from 11.5 
to 34.5 GHz, and having elevation angles from 10.7 to 57°. 

Of the 47 experiments, 18 represented two or more years of 
data (24 months or more as indicated in parentheses with the 
time interval). The attenuation values given in Table 1 for 
the 1% down to the 0.001% level of occurence were taken from 
the literature cited. 

In Fig. 7 a scatter plot of the percent deviation values 
of the simple attenuation model predictions for each of the 
47 data sets is plotted together with the mean, and standard 
deviation limits. Good agreement is obtain©' the 

model. The relatively high percent deviation: high per- 

centages of time arise from the fact that the attenuations 
are low and small deviations appear as high percentages. 

This is overcome by using an absolute deviation in dB as 
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Code 


Location 


A 28 
A73 
A 74 


Kashima, Japan 


A4K 
A49 
A40 
A 65 
A32 
A2I 
A6[ 
A 20 
A. 14 
Aiy 
A51 
A 72 
A 6ft 
\64 
A58 


A22 

BI9 

B23 
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Uluckshurg. Virginia 
(remote)]: 
f udno. Italy 
Larin. Italy 
Spinu D'Adtia. Italy 
Slough. United Kingdom 
Gometz-lurVille. France 
Leeheim. West Germany 
Kashima. Japan 
Blacksburg, Virginia 
Austin. Texas 
Waltham. Massachusetts 
Crawford Hill. New Jersey 
GrcenbeU, Maryland 
Gometz-la-Ville. France 
Slough. United Kingdom 
Martlesham Heath. United 
Kingdom 
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B27 

Palmetto. Georgia 
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Grant Park. Illinois 

* 

B32 

Clarksburg. Maryland 


B34 

Clarksburg. Maryland 
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Clarksburg. Maryland 
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i 

B31) 

Tampa, Florida 

i . 
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Austin. Texas 

.i 
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Kashima. Japan 

ih 
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4 

k ■ 
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Palmetto. Georgia 
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Palmetto. Georgia 

j ■ 
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Grant Park. Illinois 1 

j 
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Grant Park. Idioms 


C2K 

Clarksburg. Maryland 

j 
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Clarksburg. Maryland 


cut 

Clarksburg. Maryland 

| 

C2 7 
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Tampa. Florida 

! 

C2I 

Austin. Texas 
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I requency 
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(degrees) 

1 1.5 C 

46.5 

i 1 1.6 C 

107 

11.6 C 

10.7 

11.6 C 

33.0 

11.6'C 

32.0 

11.6'C 

32,0 

11. 6 C 

29.5 
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32.0 

11.6'C 

28.8 

11.7, V 

37.7 

1 1.7 C 

33.0 

11.7 0 

49.0 

11.7 C 

24.11 

11.7 C 

27.0 

1L7C 

29.0 
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33.6 

ILK C 
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ILK C 
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1 9 .04 V 

18.5 
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38.6 

19.04 V 
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19.04 V 

49,5 

19314 V 

273 

19.04 V 

41 8 

19.04 V 

21.0 

1 9.04 V 

4L<1 

19.04 V 

43.5 

19.04 V 

38.2 

19.04 V 

57.0 

19.04 V 

52.0 

19.50 C 

48.0 

28 56 V 

45.0* 

28.56 V 

38.6 

28.56 V 

29.9 

28.56 V 

49.5 

28.56 V 

27.3 

28.56 V 

41.8 

28.56 V 

21.0 

28.56 V 

41.0 
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43.5 

28.56 V 

38,2 
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57.0 

28.56 V 

52ti 
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4.15* 

34 5 C 
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Figure 7. Scatter plot of relative deviations of the simple 
attenuation model from the data of Table 1. 
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presented in Table 2. Another measure of the quality of fit 
for a model is that recently proposed by the CCIR [1981a], 

In this method, data from many radio link experiments are 
tabulated at fixed probability levels, such as 1.0, 0.1, 
0.01, and 0.001% of the year. A test variable is computed 
from the logarithm of the ratio of predicted to measured 
attenuation. To suppress measurement inaccuracies, the test 
variable is set to zero if the measured and predicted values 
differ by less than 1 dB. The figure of merit, D, is com- 
puted for each probability level from the mean and standard 
deviation of the test variables. According to the CCIR, the 
best prediction method produces the smallest D values. This 
evaluation method represents an important first step toward 
developing a standard model evaluation method. The results 
of deviation in dB and the D values are shown in Table 2 for 
the simple attenuation model (SAM) . 

Ari important consideration in the evaluation of an atten- 
uation model is whether or not the model offers an improve- 
ment over existing models. With this in mind, the perfor- 
mance of the proposed model was compared to the global model 
[Crane, 1980] and the CCIR [1981b] model. These models were 
chosen because they represent different approaches to the 
modeling problem. The global model is a rain profile model 
based on rain gauge measurements . The recently introduced 
CCIR model is an empirical model derived from terrestrial 
and slant-path attenuation data. There is a "maritime" and 
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TABLH2. Summary of Deviations Between Measured Attenuation From Data Sets of Table I and Model Predictions 
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a more complicated " continental" version of the CCIR model. 
We use the maritime CCIR model for all calculations since 
its results are superior to those found when including the 
different procedure for continental regions. The CCIR model 
for percentages of time from 0.001 to 0.1% are found in 
[CCIR, 1981b] and for 0.001 to 1% [CCIR, 1981d]. 

Attenuation values predicted using the global and CCIR 
models are given in Table 2 for the 47 experiments. The 
specific attenuation coefficients used were obtained from 
the sources recommended by each model. The CCIR rain cli- 
mate regions were used to determine rainfall statistics for 
each model to eliminate the effects of variations in rain- 
fall distribution models. Results are shown for the 1, 0.1, 
0.01, and 0.001 percentages of time, blit similar values 
occur for intermediate values. The comparison in Table 2 
indicates that on the whole all three models provide good 
agreement to data. At high percentages of time (1 and 0.1%) 
the CCIR model yields the lowest mean and standard devia- 
tions in dB as well as D values. For the important low per- 
centages of time (0.01 and 0.001%) the SAM and CCIR models 
give very good fits to the data. The global model is 
slightly inferior at all levels. The SAM model has the low- 
est percent standard deviation at all levels. Comparison on 
the basis of the D value should be made with caution. Much 
lower D values result from underprediction than from over- 
predictions by the same amount, especially at low attenua- 
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tions (high percentages of time). Thus, the D values asso- 
ciated with underpredictions at high percentages of time in 
Table 2 are disproportionately high. 


3.5 CONCLUSIONS 

In this chapter the simple attenuation model (SAM) was 
introduced for use in estimating rain- induced attenuation 
along an earth-space path. The model is both conceptually 
and computationally simple. It is a rain profile model that 
uses an effective spatial rain distribution which is uniform 
for low rain rates and has an exponential shaped horizontal 
rain profile for high rain rates. The spatial distribution 
function was derived from direct observations of rain using 
rain gauge data and verified indirectly by comparison to 
slant-path attenuation data from many experiments. Model 
estimates of attenuation as a function of point rainfall 
rate are easily computed with SAM using the physical parame- 
ters of the earth station location (elevation angle, lati- 
tude, and altitude) and the frequency of operation. To pro- 
duce attenuation exceedance estimates, the model is combined 
with a model of the point rainfall intensity distribution 
for the geographic region of interest . 

Attenuation data for various percentages of time were 
presented for 47 experiments throughout the world. See 
Table 1. Comparisons were made to this data base with pred- 
icted values from the SAM, global, and CCIR (maritime) 
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models using CCIR rain climate regions rainfall statistics. 
See Table 2. The SAii model performed well in the important 
region of low percentages of time (0.01 and 0.001%) and the 
lowest percent standard deviation at all percent time values 
tested. Furthermore, the SAM model is easy to use and is 
modular in construction. It is basically an attenuation 
versus point rainfall rate model that is coupled with rain „ 
rate exceedance to produce an attenuation exceedance predic- 
tion. This allows for separate inclusion of rain rate sta- 
tistics that affect the accuracy of attenuation exceedance 
prediction. 
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Chapter 4 

SINGLE PARTICLE SCATTERING COMPUTATIONS 

In order to study radio wave or light wave propagation 
through an ensemble of scatterers, the scattering properties 
of individual scatterers must first be determined. Knowing 
the scattering properties of the individual scatterers, the 
effects of the ensemble of the scatterers on the propagating 
wave can be studied. 

In this chapter, we will present the general formulation 
for a plane wave scattered by a single scatterer. Since the 
exact solution to this problem exists only when the scat- 
terer is a sphere (Mie-solution) , approximations to the gen- 
eral formulation will be given. After establishing calcula- 
tion procedures for scattering by a single scatterer, 
scattering by an ensemble of scatterers will be discussed in 
Chapter 5. 

4.1 THE SINGLE-SCATTERBR PROBLEM AND SOLUTION METHODS 

Let us consider a scatterer enclosed in volume v' have 
permittivity and permeability u Q . The medium in which the 
scatterer is embedded is vacuum, with parameters Eq, u q 
( see Fig. 1). E 0 is the incident plane wave upon the 
scatterer. The total electric field E(r) obeys the Fredholm 
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integral equation [1]. 
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Sc?) = E o + (vv- + k 0 2 ) f (c - £ 0 ) EC?-) Hr, ?’) av- 
where 


ip(r, r') = exp{-jk Q |r - r r |}/|r - r f 


and k 0 is the free space wave number 


k 


2 



u 


0 


£ o 


(2a) 


(2b) 


Equation (1) has an exact solution only when the scatterer 
is a sphere (Mie-Theory) . A detailed derivation for the 
scattered field by a sphere is presented in Ref. [23. For 
arbitrary scatterers, different approximations exist to 
equation (1). In most practical situations, and especially 
in propagation through precipitation, we are interested in 
the scattered field in the far zone. At a large distance r 
from the scatterer the scalar Green's function ^(r, r') may 
be approximated as 


K?, r') ~ exp{-jk Q Cr - r*r')}/r 


(3a) 
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where 
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r - 



(3b) 


and 



r 


Under this approximation the total field in the far zone may 
be written as 



(I - it) e 


r 



* 


(4a) 


whr*"- I is the unit dyadic. The field scattered by the 


scatterer in the far zone is given by 


2 -JV 


lCr) ■ 0 


k e 


4ite r 


f (0 - sj E(f’) • (I - rr) e 


ik r • r ! 
J o 


dV 1 


(4b) 


Depending on the scatterer size relative to wavelength 
and its dielectric constant, different approximations can be 
made to equation (4) . In the following, three important 
approximations will be discussed: Rayleigh scattering, Ray- 

leigh-Gans scattering and the WKB approximation. 
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a) Rayleigh Scattering [3]: 

When the dimensions of the scatterer are very small rela- 
tive to wavelength, in other words when k |r # I << 1, the 

0 

exponential within the integrand of (4) may be replaced by 
unity. Then, 


2 -JV 

k. e 0 


i(r) = ° 4lie r {p - (p • r)r> 


(5a) 


where 


P = e ojv' - 13 Efr') dV’ . (5b) 

The scatterer in this case radiates as an electric dipole of 
moment £>. 

Equation (5) holds under two assumptions £4]. Let l be 
the maximum dimension of the scatterer. Then, k 0 £ must be 
much less than unity (k Q Z << 1), and |k e r 1 | << l. 

The first assumption, i.e., k Q £ << 1, justifies the deri- 
vation of (5) and also indicates that the scatterer may be 
regarded as placed in a uniform external field. The assump- 
tion that |k £ r i | << 1 implies that the field inside the 
scatterer follows the external field instantaneously, so 
that the phase changes are of no consequence. 

The problem is thus reduced to a static one. We have to 
determine the internal field of the scatterer induced by a 
uniform electrostatic external field. By using (5b), the 
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dipole moment of the scatterer can be calculated 
internal field. In Ref. £5], the dipole moments 
oblate and prolate spheroid are calculated under 
leigh scattering assumptions. These results are 
calculate the scattering matrix of an ice-needle 
spheroid with eccentricity equal to one) , and an 
(oblate spheroid with unity eccentricity). 


from the 
of an 
the Ray- 
used to 
(prolate 
ice plate 


b) Rayleigh-Gans Scattering: 

The index of refraction of the scatterer is given by n 3 /T" 

3 

where e r is the relative permittivity of the scat- 
terer. In Rayleigh scattering, no specific assumptions 
about n were made. When n is approximately equal to unity, 
the scatterer is called diaphanaous [4]. For a diaphanous 
scatterer with In 2 - ll << 1, the Rayleigh-Gans [6] 
or Born £*7] approximation holds. 

In the Rayleigh-Gans approximation the field E(r) inside 
the scatterer is approximated with the incident field E 0 . 

Under this assumption, the far-aone scattered field § s (r) 


from (4) becomes 

2 ' j V 

k/ e 0 




4ttg r 
o 


L t. - ... tf, - 


dV' 


( 6 ) 


For a homogeneous diaphanous sphere (see Fig. 2) the scat- 
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tered field takes the form 

7 /• ik '' r 

E s (r) - k Q r - o- o 


. 2 | (e - e ) {E - (E . r)r} 

o 4ire r o o o J 
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*2tt e IT f> 3. <p (7a) 

d<j>' sin e* de* r' exp{2j k Q sin(e/2) cos e T } dr* 

o J 0 0 





. » sin(2 k sin(e/2)a) 

r)r> [ — 2 k f ' sin(9/ ' 2l ' 
o 


a cos(2 k o sin(e/2)a) 
(2 k 0 sin(8/2)) 2 


(7b) 


where a is the radius of the sphere and 0 is the angle bet- 
ween r and the z-axis, assuming that the origin is the cen- 
ter of the sphere. For more complicated shapes of scatter- 
ers the volume integral in (7a) cannot be calculated 
analytically and numerical solutions must be used. Exact 
evaluation of the volume integral can be done for ellipsoids 
[4], 

c) High Frequency Scattering; The WKD Method: 

The WKB approximation is applicable to cases where the 


Rayleigh or the Rayleigh-Gans approximations cannot be 
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applied. Specifically the WKB approximation holds when 


|n^ -l| k Q £ >> 1 and |n^ - l| < 1 


( 8 ) 


In the WKB approximation and the field E(r) inside the 
volume V' is approximated by a plane wave propagating in the 
same direction, as the incident field, with propagation con- 
stant equal to 


k2 ■ “ 2 w 0 e o = Cn V 


(9) 


Under these assumptions equation (4) takes the form 

2 0 * e o> * -2t k 0 z x + k 0 nCz ’ * z l )] 

0 


2 -JV f 2(e - e o } * 

S(r) - k 2 |— r- I n + 1 E o 
O ^Tre^r Jyr 


tlO) 


xe^o 1 


r' 


dV' 


The incident plane wave is assumed to be propagating in the 
s-direction. is the value of z' associated with r # for 

points lying on the surface of the scatterer as shown in 

2 

Fig. 3. The factor — — T is the normal incident transmission 

n + 1 

coefficient from a vacuum to the medium of the scatterer. 
The WKB approximation is part of a general mathematical 
method developed by Wentzel Kramers, and Brillouin [8]. 
According to the WKB methou the field inside the scatterer 
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is expanded in the series 
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ICr) = e jk ° f ^ {S o + + ~7 * 

0 K o 


(U) 


Using the above equation in conjunction with (2) and col- 
lecting together like powers of k Q , a set of equations can 
be obtair^d for f(r), E 0 , E 2 , ... * etc.. For 

high frequencies, k Q (n 2 - i)a >> l. in this case t(r) 
can be approximated by 


- jk o £ ^ 

E(r) - e 


( 12 ) 


and equation (10) is obtained. 

The three approximations discussed so far are valid only 
in specific cases when certain assumptions hold. However, 
in many scattering problems it is necessary to calculate the 
scattered field when none of the aforementioned approxima- 
tions hold. For these specific problems exact analytical 
results do not exist. Most of the work in this area is 
based on numerical techniques and the use of fast computa- 
tional machines. In the next few paragraphs we will summar- 
ize the most commonly used methods to calculate the scat- 
tered field by a dielectric body with complex permittivity. 
The majority of these methods have been used to calculate 
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the scattering coefficients of a spheroidal or ellipsoidal 
scatterer . 

i) Point-Matching Technique: 

In the point-matching solution the incident field, the 
scattered field, and the field inside the scatterer are 
expanded in terms of spherical wave functions. The nfinite 
modal summations of the expansions are truncated at some 
model index (M, N), and boundary conditions are applied to 
the same number of points as the number of unknown expansion 
coefficients. The coefficients are then calculated by 
inverting a square matrix that is formulated by the 
boundary conditions at M*N points. 

The method described above has been used by Oguchi [9]. 
Morrison and Gross [10] also used the point-matching techni- 
que with a least squares fit. In the latter case the number 
of boundary points is larger than the number of unknown 
coefficients and the fields are matched at these points in 
the sense of least squares. The least squares fit technique 
coverges much faster than the one used by Oguchi . 

ii) Spheroidal Wave Function Expansions: 

The scattering properties of a spheroidal scatterer are 
obtained by expanding the fields in spheroidal vector wave 
functions and truncating the infinite summation of the 
expansion to a finite sum. The unknown coefficients of the 
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expansion functions are evaluated by applying boundary con- 
ditions. Usually, spheroidal wave function expansion meth- 
ods require the boundary conditions to be applied to fewer 
points than the point matching technique, especially when 
the scatterer is a spheroid. This method has been applied 
by Oguchi [11]. 

iii) Waterman's T Matrix Formulation (Extended Boundary 
Condition) : 

This technique has been used extensively for scattering 
by perfectly conducting bodies. Recently, the technique has 
been applied to dielectrics. The scattered field is 
expressed in terms of electric and magnetic surface cur- 
rents. These currents are expanded in terms of 

[12] spherical harmonics. Boundary conditions are applied 
on an inscribed sphere inside the scatterer. By truncating 
the infinite expansion of the currents, the expansion coef- 
ficients can be obtained by matrix inversion. 

iv) Fredholm Integral Equation Method: 

This method was introduced by Holt, Uzunoglou and Evans 

[13] . Starting with the integral equation for scattering of 
(1), it is shown that the Fourier transform of the field 
inside the scatterer is the solution of two coupled integral 
equations. The integrations are reduced by numerical qua- 
drature methods to matrix equations, whose solution can be 


- 60 - 


easily obtained. It is important to notice that the scat- 
tering amplitude obtained by this method satisfies the 
Schwinger variational principle, and thus the method is sta 
ble. 

Of these numerical methods the most favorable are Water- 
man's method and the Fredholm Integral Equation method, 
since both converge rapidly over a wide rcinge of scatterer 
sizes. It should be noted though that the latter needs a 
large amount of computer storage when the shape of the scat 
terer is complicated. 


4.2 THE FREDHOLM INTEGRAL EQUATION METHOD APPLIED TO 
RAINDROPS 

An extensive derivation of the method is described in 
[13]. The details of the derivation have been verified, but 
in this section we shall present a summary useful for per- 
forming calculations. Example calculations are performed 
for rain. The size of raindrops is of the order of the wav- 
elength in the microwave region and Rayleigh or Rayleigh- 
Gans approximations do not hold. 

The field scattered by a dielectric body of relative die- 
lectric constant e r (r) and volume V obeys the equation 


y(r') 


(13a) 


where 
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G (r, r 1 ) = 


1 + 


l - 7 7 
2 - - 


exp(i k Q r - r f | ) 
4 it I ? - ?' I 


(13b) 


y(r) = k 2 [e (r) - l] (13c) 

o r 

k Q is the free space propagation constant, 1, is the unit 
dyadic, kj is the direction of propagation of the incident 
wave and 


J, = 1 - k. k. 
=i = ii 


(14) 


In the far- field of the scatterer (r ® ) the scattered 
field in the direction kg is given in terms of the scat- 
tering tensor f(k s , k ± ) as (see Fig. 4) 


exp(i k r) 


E (r) - 

3 


f= ( W 


(15a) 
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Figure 4. Geometry for electromagnetic scattering 
from a single spheroidal particle. 



where 
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1(5 » 5. ) 

— s 1 


1 

4ir 



exp(- k s 


5) y(5) E <?)dr\ 


(15b 


Let C(5> be the Fourier transform of E(r) with respect to r 
Then f(k , k.) can be expressed in terms of £,(k) as 

3 1 


£<V V ■ 4^ 4 ■ / d k 2 v ^s> V £ ( V 


(16a; 


where 


V(5 X * 5 2 ) y(r) exp£-i(5^ - 5 2 ) • rjdr 


(16b; 


C(5) obeys the integral equation 





• C (k 2 ) = 1. V(k r 1) 


(17 


where 


£,) = IwCk, 5„) - ICk-, k„) 


(18a 
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and 


m v k 2 ) = 


lim 


/- 


2 j 

p d p 


**«. j 2 2 
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I-PP 


V(^,p) V(p,$ 2 ) (18c) 


Evaluating the integrals o£ (15b) and (17) by numerical qua- 
drature we get: 


£<V fy = h a ( V v <>V fy U9 ' 

where 

S w & Uk., t^) • £(k^) - ^ V(kj, j=l,...,n (20) 

^“1 

Equation (20) can be used to solve for ) w a by matrix 

inversion? and then placing these values into (19), the 
scattering tensor of the dielectric scatterer can be calcu- 
lated. 

A spheroidal dielectric scatterer obeys the equation 



(21a) 


i 


1 
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with 
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b = a ( 

-y — |> 

Assuming that e r (r) is constant and equal to e r , V(kj_, 
k2> for the spheroid becomes 


VCk^, Ic^) - 47rabc y j ^ CX) /X 


with 


X = 



and 


jn(x) 


is 


k. = k. (a sin 0. cos <b , , b sin 0. sin <p . . c cos 0,) 

ii i i i T i i 

the spherical Bessel function given by 


j ( x ) = (JL ) 1 / 2 r 

V x; ^2x ; * (n + 1/2) 

Also the tensor j£{ki, Tc 2 ) for spheroidal shaped scatter- 
ers takes the matrix form: 


K(it, £ ) -el t ) - &(S , t ) 



where I is a 3x3 unit matrix and 
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(30) 


X i - cx[a 2 + (c 2 - a 2 )x 2 ] ^ 2 i=l ,2 ( 31 ) 


A»t < x > 
1/2 
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and 


Sst 00 


2 

I^s.ti-y I fi (s,e) 
-y 2 I 5 (s,t) 

-xy I 3 (s,t) 

\ 
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-y 2 I 5 (s,t) 

-xy 

I 3 (s,t) 

X 1 (s,t)-y 2 I 4 (s,t) 

-xy 

I 2 (s,t) 

-xy I 2 (s,t) 

2 

y 

I^(s, t) 


J (32) 


Notice that the matrix Q s t(x) is symmetric and from (26) 
it is obvious that K(kj., k 2 ) is symmetric. The elements 
of the Qg t (x) matrix are given below. For s+t even 


I, (s,t) = G(4>. » 


i, = r = o 


i. 

f \ 

J. 4. 

4. O 

/ 


f ' 

I 4 (s,t) 
< I g (s,t)’ 

> = G(<t> r 4> 2 ) ( 

sin 2 (<}>) 
sin <t> cos <p 

> + * { sc < 

cos (s(Ji^ 2 +4 1 ^) 
-sin(s0 12 + ) 

I 6 (s,t) 


2 T 
cos <p 


-cos(s4i 1 7 + 24 ) 

l ) 

< J 


K ' 
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and for s+t odd 


I 1 (s > n) = I 4 (s,t) = I 5 (s,t) = 1^ (s, c) - 0 
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I 2 (s,t) 


I 3 (s,t) 


GC^* 4> 2 ) 


sin <f> 


cos 4> 


(34) 


where 


* 2 > = 2 Cj|j < [cos 4> 12 ] 


(35) 


C m r (x) is the Gegenbauer polynomial 


’12 = *1 " *2 


s ^ t 


(36) 


s < c 


(37) 


In the version of (37) presented by Holt et al , the equal 
sign was missing in the inequality n > m. 
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i # 0 2 are the angles of wavevectors k 1# k 2 , respec- 

H 

Lively. See (24). In (27) p n (x) is the associated Len 
gendre function and h^fx) is given by 


h n<*> * <^ 1/2 H nil/2 <*> 


(38) 


where H n+ ]/ 2 (x) is the Hankel function. After finding 

-H *> 

K(kj # k^) where j»l, 2, ...,n£ = 1,2, ...,n and 


-v -»■ 



k. ) we can solve (20) for 


Cc(k^, C(k 2 ) ... , C(k n )]T, The scattering matrix 
■+■ ->■ 

f(k s , k^ ) can be computed then from (19). 


-*■ -*■ 


In calculating the element of the Z(k^, k 2 ) matrix 
the spherical Bessel and Neumann functions must be calcu- 
lated as well as the associated Legendre functions and 
Gegenbauer polynomials. The integration in (27) required to 
evaluate the Z(k^, k 2 ) matrix was performed using a 
seven point Clenshaw-Curtis quadrature. The infinite summa- 
tion was truncated to N terms. It was found that conver- 
gence occurred for small raindrops when N = 5 and for large 
raindrops when ^13. Note that these values of N were also 
used by Holt et al. The wave vectors k^, k 2 are com- 
posed of ( x^ = cos ) and (x 2 = cos 0 2 , $ 2 ) 

respectively. In (20) we chose the x pivots to be n x with 
x varying between -1 ^ x < 1. The n x pivots are equally 
spaced in the above interval. The pivots are equally spaced 
in the interval 0 < <f> < it. 
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A Fortran computer program {called COEFF) was written. 

In Section 6.2 the COEFF program is described and a state- 
ment listing is presented. 

Finally in this section we present the results from an 
example computation for raindrops using the COEFF program 
and compare those results to published values. Let 

A A. 

f f (0) = 8 (k Q x, k Q x) . z ( 39 ) 

f"(0) = S (k Q x, k Q x) • y (40) 

f'(0) and f" (0) are the forward scattering coefficients for 
incident field polarizations along the major and minor drop 
axes, respectively. The values f'(0), f* 1 (0) were calculated 
for a spheroidal raindrop with a = 0.125443 cm, c = 0.172355 
cm and s r 2 - n Q = 7.B84 + j 2.184 at frequency of 11 
GHz. These values were compared with one's of Holt. (See 
Table 1) For the above calculations the pivots n x , were 
both equal to 5 and N=9. The execution time was 10.5 
minutes. Execution time increases significantly with 
increasing N, n x , and n^ . Keeping n x , constant, CPU 
execution time is proportional to W2 . p 0 r larger values 
of a, N must be increased which causes a large increase in 
execution time. 
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4.3 THE SINGLE-PARTICLE SCATTERING TENSOR 
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We have seen that single-particle scattering coefficients 
can be calculated by one of a variety of analytical or num- 
erical methods. These results are then to be used in com- 
puting the scattering effects of an ensemble of particles in 
precipitation media. This is facilitated by casting the 
single-particle problem into a tensor form. In this section 
we discuss the single-particle tensor in preparation for use 
in the multiple scattering algorithm for rain in Chapter 5. 

For an incident plane wave propagating in the z-direction 

— w 

with electric field E 0 , the far zone scattered field of 
(4) has electric field components perpendicular to the 
direction of scattering and can be expressed as 


-K ,-K 

V r > 



£<r) E o 


(41) 


where f(r) is the scattering tensor of the scatterer. For 
example, in the case of Rayleigh-Gans scattering it follows 
from (6) that the scattering tensor is given by 


f(r) = 



o 


4tte 

o 



_ +jk r • r' 

e ) (I - rr) e 0 
o 


dV' 


(42) 
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Equation (41) may be written in a matrix form as 


- - 


" 1 

E 

X 

s 

-jkr 

- 5-j— f (r) 

E 

ox 

E 

y s 


E 

oy 



a m 


where £(r) is a 2x2 matrix corresponding to the scattering 


tensor 



f(?) = 





-J 


(44) 
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Chapter 5 

MULTIPLE scattering computations 
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5.1 Introduction 

Scattering problems involving more than one scatterer are 
complicated and, in general, analytical solutions do not 
exist. It is necessary to make certain simplifying assump- 
tions. Typically it is assumed that the scatterers are ran- 
domly distributed, are infinite in number, and are in the 
far zone of each other. In this paper a, general approach to 
scattering from particulate media which includes multiple 
scattering is presented. The results are applied to the 
problem of millimeter wave propagation through rain. An 
excellent review of the role of multiple scattering in 
radiowave propagation through precipitation is given by 
Olsen [1982]. 

Scattering by random distributions of scatterers was 
first studied by Rayleigh [1871], who used a single-scatter- 
ing approximation and identical, aligned scatterers. A 
detailed derivation of Rayleigh's results is available [Van 
de Hulst,1957]. The single-scattering approximation does 
not hold when the density of scatterers is large or the pro- 
pagating wave frequency is high. Improved accuracy over 
that using single-scattering is possible using first-order 
multiple scattering or complete multiple scattering. 

In this chapter equations are derived for the average vec- 
tor electric field <E(r)> in the presence of a random dis- 
tribution of scatterers using, first, lower order scattering 
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approximations and, second, the multiple scattering 
approach. In Section 5.2 the single-scattering approximation 
is used to derive expressions for the coherent electric 
field in an ensemble of scatterers with random distributions 
of location, size, shape, and orientation. These results 
are used to treat the ensemble of scatterers in which parti- 
cle scattering only occurs once but the incident wave on 
each particle has included the effects of previous parti- 
cles; this is first-order multiple scattering. In Section 3 
complete multiple scattering is considered in which the 
field incident on a particle can arise, in part, from fields 
scattered from any other particles. The Foldy-Twersky scat- 
tering procedure is used to derive an integral equation for 
the vector coherent field in an- ensemble of particles. In 
Section 4 multiple scattering results are applied to the 
rain propagation problem. Some typical calculations are 
performed. 


5.2 Lower Order Scattering 

Scattering from discrete media that does not include all 
orders of multiple scattering has been studied from various 
approaches by many investigators. An excellent review was 
prepared by. Ishimaru [1977] . In this section we discuss two 
of these which have been applied to the rain propagation 
problem. 
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Single scattering 

Under the single-scattering approximation, the electric 
field is scattered only once by the scatterers. Let us con- 
sider specifically the slab of scatterers of thickness L 
shown in Fig. 1, with a plane wave incident upon it. Each 
scatter er (say the p* 1 * 1 one) can have’ a random oosition 
r p , size as measured by the equivolumetric radius a p , 
shape as represented by the shape parameter s p , and orien- 
tation angle represented by 0 p . The field at point r, 
where the receiving antenna is located, is the incident 
field plus the field scattered by particles in the active 
region of the slab; the active region is the central few 
Fresnel zones. Summing over the contributions from W parti- 
cles in the active region, we obtain 


E (r) - E 1 + 


N 

l 

p=i 


k (x ' 
o v p 


y p 2 )/2 


I (?) 

p 


E 1 (r) 


( 1 ) 


where f D (r) is the scattering tensor of the p fch scat- 
terer and r p = |r p | . 

In order to compute the average field at point r we 
assume that the scatterers are randomly distributed in posi 
tion, size, shape, and orientation angle, and that all of 
them have the same particle distribution 



) 


( 2 ) 
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where Sp represents the random parameters of position 
rp, size ap, shape s p , and orientation angle 9 p . n(S D ) 
is the number of scatterers per unit volume for parti- 
cles in class w n . 


Equation (1) is averaged using the distribution given in 


-* _ 

<E Cr 3 > ~ I + 

_ ■’ r ' 


7(r f ) 


-jk 0 [x' z + y ,2 )/2r 


n (r 1 , ui) dr * dw * E 


Here, dr 1 - dx f dy T dz ’ is the elemental volume of space, and w 
encompasses the distribution parameters for the particle 
size, shape, and orientation, w = (a, s, 0) . 

Since major contributions arise from particles in the 
first few Fresnel zones, r* is nearly independent of x*, and 
y' in in the integral appearing in (3) , We can make then 
the substitutions 


x i = 




- v’ 

r / J 


= Z' (4) 


giving 


2 , ,00 CO 

I 1 

<£ITJ> = 1 + Jj7- dw dz, dx.. dy^ 

J n J J J ~ a? < -» 


-Cx i 2+y i 2 ^ - 
e f nCr-^j w) 


->4 -y 

• E 1 (r) . 


1 


Using the relationship 

00 „2 

e u du - /tt 

j. co 


0/r POOft 


FVSr,-*.' - 

_ ' to 

Quality 


( 6 ) 


and assuming that n(ri* w) varies only in the direction of 
propagation z, (5) gives 



(d) dz 



■*i 

E X (r) 


- [I - jlc] • E i (r) (7) 

where 

k = y ~ j T n(z^, «) dZjdw . (.8) 

If the elements of the tensor k are small, the quantity I - 
jk may be approximated by an exponential. This step is 
further motivated by our knowledge that the exponential form 
includes multiple scattering effects as it will be shown in 
Section 5.3. Within the framework of this approximat >n (7) 
can be written in the matrix form as 


<E(r)> = e'jjiEMr) (9) 

Here, <E(r)> is the column matrix corresponding to the vec- 
tor <E(r)>. Similarly, E-^r) is the column matrix associ- 

-*• ■ -9“ 

ated with E 1 (r) and _k the matrix corresponding to the ten- 
sor k. 
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The coherent field scattered by a slab of scatterers can, 
in general, be written as 


<E s (r) > = D * E 1 (r ) (10) 

where D is the scattering tensor of the slab and depends on 
the scattering properties of the individual scatterers and 
their random distribution in position, size, shape, and 
orientation angle. Specifically, for the single scattering 
case 


0=1 -jjj 1 ?(to) n(z 1 , w) dZjdw . (11) 

O * Z £ , U) 

In the derivation of this result the assumption was made 
that the plane wave is scattered only once by the particles. 
This assumption is unrealistic, since the wave may be scat- 
tered several times by the particles before reaching the 

+ 

point R at position r. When the scattering coefficients of 
the particles and the particle density are small, contribu- 
tions to the scattered wave by second, third, and higher 
order scattering may be ignored. There are cases, however, 
where the scattering coefficients of the scatterer are large 
(e.g., for raindrops above 30 GHz), and the multiple scat- 
tering contributions cannot be ignored. Starting with the 
next section we will take into consideration multiple-scat- 
tering contributions. 
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First-order multiple scattering by a slab of scatterers 

Under the first-order multiple scattering approximation, 
the slab of particles of length L is divided into n thin 
subslab's. In each individual subslab the single approxima- 
tion is assumed to hold. The coherent field at point R is 
then given by 


<E(r)> = D 1 ■ D 2 • ... D n • Ei(r) (12) 

where is the scattering tensor of the i fc ^ subslab and 
is defined by {11} . This first-order multiple scattering 
approach has been used by a number of researchers in differ- 
ent areas of wave propagation. In radio wave propagation 
through precipitation, for example, Persinger et al. [1980] 
have used a computerized code to evaluate (12) in order to 
calculate attenuation and isolation of a wave propagating 
through a medium consisting of ice-particles and rain. 


5,3 ?iultiple Scattering 

Single scattering holds only when the scattering coeffi- 
cients and the density of the particles are small. First- 
order multiple scattering holds for any forward scattering 
propagation situation provided that the thickness of each 
subslab is kept small. As shown in (12) , the coherent field 
in this case contains the dot products of a large number of 
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tensors and its evaluation can be performed only by fast 
digital machines. 

In order to include both forward and backward multiple 
scattering, we will use an approach first introduced by 
Poldy [1945] and developed further by Lax [1951] and Twerskv 
[1962] . In our formulation we will assume that the medium 
is anisotropic, and we will derive the Dyson equation for 
the coherent vector electric field in terms of the scatter- 
ing tensor of the individual scatterers and their distribu- 
tion in location, size, shape and orientation. 

Derivation of the general multiple scattering integral equation 

Let N scatterers be randomly distributed in space, and 
have random size, shape, and orientation. The medium bet-' 
ween any two scatterers is free space. (See Pig. 2) The 
vector electric field E(r) between the scatterers satisfies 
the Helmholtz equation 

(V 2 + k Q 2 ) E (r) = 0 , (13) 

where k Q is the free space propagation constant. 

The incident field E a at the point r a is the sum of 

+ _ 

the incident electric field Ei a in the absence of the 
scatterers and the field scattered by the N scatterers, spe- 
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the incident field in the ubsen 


He crnttf'rf'd -F-rnm all n1*hpr 


na 







cifically, we have 
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N 

E a = Ei a + l Gj a , (14) 

j*l 

+ -> 
where Gj a is the field at position r that is scattered 

by the j th scatterer. The latter is a function of the 

scattering properties of the j th scatterer at position 

?j and the location r. We can define a tensor gj^ 

such that 


Gj a = gj a * 


The field incident on the jth particle is 


= eJ 


N 


• J + l G.- 
1 k=l,k^j k 


Substituting (16) into (15) we obtain 


“►■j a “*■ i -*■ 4 

G J • g/ . [E. J + l G J] . 

] 1 k=l,k?j K 


(15) 


(16) 


(17) 


Using (17) and (14) , we nave 



N 

I 



E i 


j 


N 

+ I 
j-1 


N 

l 

k=l,k^j 




E k 


. ( 18 ) 


■ f 


^'OINAL P Aq , : ... 

By iterating the above process it follows that *■ P °°« OtlAUiv 


, . IN _W * IN iN _ • .+, 1, 

E a = E, a + I f, a • E, 3 + l l g, a • g k J • E^ k 

1 j = l 3 1 j=l k»l,k?fj 3 k 1 


N N N 

Z r Z 

j =1 k-l,k?*j Jt = l,£^k ** 


— k c £ , 

• • E. + . . . 

s £ 1 


The triple summation in (19) may be written as 


S ? ? -a - j -k * A 

I I I Si • g k • St • E. 

j=i k=i,Mj i-i,m 3 k 1 


N N N _ 3 _ i _ lr - » 

* l l l gi • g k J • g^ k • E a 

j = l k=l,ki<j i-l,i«k,i^j 3 K 


N N _ , _ k * , 

* l I gj 3 • g k J ■ gj k ' E. 3 

j=l k=l 3 K 3 1 


( 20 ) 


Assuming that backscatter ing is smaller than forward 
scattering, the second summation in (20) may be ignored, and 
(19) may be written as 


E = E 


„ in „ • N W _ 0 -> 

. a + l g. a - E. 3 + y j g. a * g a * E. 

1 jii §3 1 jh k-itwj §3 Sk 


« W « „ a _ i _ k -► 

+ I I I * g k J * g £ k • V 

j =1 k^l^j £«l,£j*k t £j*j J 


Another way to state the significance of the approximation 
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involved in this equation is that we have ignored the 
following: triple scattering between two scatterers, fourth 

order scattering between three scatterers, and so on. The 
above procedure was first introduced by Twersky [1962] in 
order to obtain a closed form equation for the coherent 
field. 


We next make the assumption, as in Section 2, that the 
scatterers have the same position, size, shape, and orienta- 
tion distribution, and that there are no correlations bet- 
ween scatterers. The particle distribution may be defined 
as in (2). However, now N is the total number of particles in 
the entire volume. Taking the ensemble average of (21), we obtain an 
expression for the coherent field at point r: 


<E a > 



N 

- i 
j=i 





N j 


N N 

- l l 

j=l k=l,k?ij 



-► k n(u.)n(u.) 

■ P _ mu . ^ — 


* * N N N 

E i“ ,u ^2 — — dw. du k + l l l 


N 


j-1 k-l,k^j 




-a -k t t n( Mj )n( Mk )n(M t 3 

O' . • CT • h . -..r' J , 

N* 


g k" • V • V — J — zr — ~ dS j d “k d “i + 


(22) 


or 
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<E a > 


“i # * gj 3 ’ E.J a (a.) dZ j * Sji | 




m 3 


— a — i c 

Si ’ S k J • E i 


nC*',)n((5 fc )dS, dS k + ^ L + * S 


“j ,w k 

— a — j — k t i 
’■ ’ g k • *» * E i 


nCwj )nC^)n(aj Jl ) daK du^ dtu £ + .... * 


(23) 


By letting N tend to infinity, the ratios — , etc,, 
tend to unity, and the infinite summation in (23) can be 
represented by the integral equation 

<E a > = E- a + f g. a * <E^> n(w.) dai. . (24) 
x JJ, ^ 3 3 

This is the Dyson equation for the coherent field in an 
ensemble of randomly distributed scatterers, and it is often 
called the Foldy-Lax-Twersky equation. It has been derived 
on the basis of the Twer sky procedure. The same equation 
may also be obtained using Lax T s polycrystal ine approxima- 
tion [Lax, 1951] . 


The Dyson equation, (24) , is more general than the equa- 
tions for the coherent field derived by the single scattering 
or first-order multiple scattering approximations? it 
includes backscattering. Evaluation of (24) is very diffi- 
cult and depends on the complexity of the operator gj a . 

If the medium is tenuous (the average distance d between any 
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two arbitrary scatterers much greater than the free space 
wavelength, i.e., k Q d >> 1), then each scatter er is in the 
far zone of all other scatterers, and gj a may be 
replaced by 



C25) 


where f is the scattering tensor for the scatterer. 


Evaluation of nriiltiple scattering for plane wave incidence 


Most practical situations {such as in communications 
applications) are well approximated by a plane wave incident 
normally on a slab of scatterers as described in the previ- 
ous section. Let a plane wave be propagating in the z-di- 
rection with electric field 

-*■ + -J*o z 

E. = E.e ° 

1 o 


The average field <E> inside the slab obeys the vector Fol- 
dy-Lax-Twersky integral equation of {24) . The medium of the 
scatterers is assumed to be tenuous, so that the scattering 
operator gj a can tahe the form shown in {25) , Substi- 
tuting (25) into (24) and evaluating at the observation 
point a with position vector r, 


<E(r)> 




z 


L fjCz, r-fj)<S- l > nCu 3 ) d£. 

J oi. ■' J ] r - r j | J J 


+ 


(26 


ORIGINAL 

where |r-?j| - [(x-Xj) 2 + (y-yj) 2 + (z-zj) 2 ] 1/2 0F POOR quality 

and tuj * (rj aj, Sj, 9j). The medium is assumed to 
vary only in the z-direction and, therefore, the fields will 
also. Evaluation of (26) then proceeds by noting that the 
medium is assumed to be homogeneous in the x,y-plane and 
employing the method of stationary phase to reduce the 
integral giving 


+ *> -Jk.* 2 it i 

<E(z)> = E e 0 + 4^1 

° k o 


_ •» -jk Cz-z>) 

7 <E(z’)> e 0 n(u>) doi (27) 


CD 


where ui = (z f , a, s, e) and f is the single-particle forward 
scattering tensor. The solution to this integral equation 
in matrix form is 


<E(z) 


rii 


E e 
o 


- jk o 2 


(Z8a) 


where 

— ” f — n(<*0 dw C 2 8b ) 

o ai 

& similar solution has been reported [Oguchi, 1981] in the 
case where the medium was uniform in the z-direction. 


Note that .'28) is the same as 1 9 ) • Thus, the exponential 
approximation to the derived single-scattering formula in 
(7) is justified on the basis that it really includes multi- 
ple scattering effects. 


The equations derived in this section for the coherent 
vector electric field can be used to study radio wave propa- 
gation through precipitation media which vary in composition 
along the direction of propagation. This is explored in the 
next section. 


5 . 4 ^ipli.cations to Rain Media 

The results of the derivation in the previous section can 
be applied to precipitation media. Here we consider rain. 
The formulation has the following features: it includes 

multiple scattering effects? it allows for direct inclusion 
of distributions of raindrop sizes, shapes, and canting 
angles; and it allows for a nonunifornr rain rate along the 
propagation path. Most models used to make actual calcula- 
tions for rain do not have such flexibility. Quite fre- 
quently the rain is assumed to be of uniform rate over some 
path length, the canting angle distribution is omitted or is 
included by modifying a constant canting angle result, 
and/or all drop shapes are alike. 

Persinger et al. 11980] presented a first-order multiple 
scattering model that included a nonuniform spatial rain 
rate variation as well as raindrop canting angle and shape 
distributions. In this section we use the ideas of Per- 
singer et al. together with the multiple scattering model 
to make calculations for rain media. 
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The equivol ume trie radius of the raindrops varies from 
0.1 to 4.0 mm. At millimeter wave frequencies the raindrop 
radius is comparable with the free space wavelength and the 
Rayleigh scattering approximation is not adequate for the 
calculation of the raindrop scattering matrix* It is neces- 
sary then to use one of the available numerical methods 
[Oguchi, 1981] to calculate the raindrop scattering coeffi- 
cients. In our computations we will use values of scatter- 
ing coefficients published by Uzunoglu et al. [1977] . 

These results have been computed using the Fredholm integral 
equation method. We will assume that the raindrops are 
either spherical or oblate spheroids. 

Let e be the canting angle of the oblate raindrop.' The 
single raindrop scattering matrix is given by 

f v (a)+f H (a)+cos 2e [f H (a)-f v .(a)] [fy(a) -f H (a) ]sin 2 e 

[fy(a)-f H (a)]sin 2e f v (a)+f H (a)+cos 20 [fy(a)-f H (a) ] 

In this expression r a is the equi volumetric radius of the 
oblate raindrop, and f v (a) , fg(a), are the scattering 
coefficients for the incident electric field aligned with 
the minor and major axes, respectively. The quantities 
f v (a) , f H (a) are expressed in terms of powers of the 
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equi volume trie radius of the raindrop? specifically, 

5 i 

£Vfa) - l a, V a 1 £30) 

H 1=0 1 H 

Values of ct.V are obtained by curve fitting published tabu- 
X H 

lar results in terms of size, at a specific elevation angle. 
Rain medium computations 

The rain medium is represented by a slab of raindrops of 
extent L. With plane, wave incidence the coherent field pro- 
pagating inside the rain medium obeys (28), where the single 
drop scattering matrix f is given by (29) . In nCz 1 , oj) 
n(o0 of (28) , z* is the distance along the direction of 
propagation into the rain slab, a is the equivolumetric 
radius for the individual raindrops, s is the shape of the 
raindrop (spherical or oblate speroidal) , and 0 is the rain- 
drop canting angle. In the distribution density function nCw) 
of (28b) 0 and s are assumed to be statistically independent; 
i.e, . 


n(oj) * n 0 (z», a) p x (9) p 2 (s) . (31 ) 

P 2 (s) is a discrete shape distribution density function in 
which P 0 is the fraction of drops that are oblate spheroi- 
dal and the remaining fraction 1 - F 0 are spherical. The 
canting angle distribution is Gaussian with mean <e> and 


standard deviation a 

PiO) 


9* 

1 


- Ce -< 0>) 2 /2 a Q 2 


✓27 
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0 


For n 0 (z, a) we use the drop size distribution of Marshall 
and Palmer [1948] with rain rate as a function of positions 


n 0 Oi a) = 16,000 e “8-2LR(z) °' ZI ]a 


-3 -1 

in mm 


(33) 


The rain rate spatial variation used in the examples to fol- 
low is piecewise uniform as introduced by Persinger et al . 
[1980] , where 



0 < 2 ^ 0.2L 
0*2 L < 2 *£ L 


C34J 


in which R 0 is the point rainfall rate at the receiving 
station (z - L) . 


Evaluation of k in (2Sb ) can proceed after substituting 
in the above stated information on the distributions: 

£ - f* dz ’ f \ d9 Ul-F 0 3 f SPfi Ca)+F 0 £ 0BL Ca,e)}n o C 2 ',aj ri (e) 
The integral over the angle 0 may be performed by first 


C 
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evaluating the expressions 


J < 0>+ct 

do sin 2e (6) 

<6>-a 


<cos 20> 



<e>+ot 

<0>-a 


de cos 20 p^(9) 
(3 


where a is the maximum variation of the canting angle. It 
can he shown that iC36) is approximated by 


<sin 20 > ■ exp(-2cr Q 2 ) sin 2<e> 


(37 ) 


and, similarly, 

<cos 20> » exp(-2cr e 2 } cos 2<9> . (33) 

With these approximations the average scattering matrix k 
becomes 

fc = § - dz' f da{Cl-F 0 Jf SPH La) + F a <f 0BL Ca]>>n o Cz,a) C33) 

where <fOBL( a )> the average of fp SI, (a, q) with res- 
pect to 0, and is obtained from f(a, 0) by replacing sin 20, 
cos 20 with <sin 20> and <cos 29>, respectively. 

OBL 

The scattering coefficients f s ^ H (a), £V (a) are 

. H 

expressed, as in (30) , in terms of powers of a up to the 
fifth order. Then the final evaluation of (39-) follows 



after performing the indefinite integral 


I 



j 


da f (a) n o Cz f , 


a) 
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- I c_ f L dz 1 [ da a n e' YCz)a C40) 

n=»0 n J o J 

where y (z # ) = -8.2 an( j Cn are the expansion 

coefficients as in (30). Integrating over a and substitut- 
ing in the rain rate variation of (34) and performing the 
final integration gives 


i 





i 


h 

c 

I; 

0 


2 5 n I 

- I *i I c n I 

1=1 1 n=0 n m=l y. A 


r _ n — in n 

- 1 * 3 C-D® + ~ 

Y i 


(41 ) 


where = 0.2L, S>2 m 0.8L, and =* Y(&£>. 


The elements of the matrix h have been evaluated explic- 
ity with (35) and (4tL) . The depolarization matrix e*i!s of 
(28a ) is then evaluated using eigenvalue-eigenvector meth- 
ods. After the depolarization matrix is evaluated for the 
rain medium, attenuation and channel isolation on a dual-po- 
larized system can be calculated for a wave passing through 
the. rain slab [Cox, 1981]. 


Comparison of rm J.t±ple and first-order multiple scattering resists 

The multiple scattering model as presented in Section 5.3 
has been computer programmed. Program inputs include fre- 
quency, elevation angle, mean and standard deviation of 



\ 

) 

i 

i 


i 

i 

i 
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canting angles , rain rate at the receiver, fraction of 
oblate raindrops, rain medium length, and incident wave and 
receive antenna polarization parameters. 

Calculations have been performed to test the model and to 
examine the frequency dependence of multiple scattering 
effects. To facilitate comparison, as frequency is varied 
over 11, 20, and 30 GHz the rain medium parameters are fixed 
with a 5 km length, a uniform rain rate along the entire 
path, and a shape mixture of 60% oblate and 40% spherical 
raindrops. The oblate raindrops have a Gaussian canting 
angle distribution with a mean of 0 degrees and a standard 
deviation of 12 degrees. The incident wave polarization is 
linear with a 45 degree tilt angle. The elevation angle is 
45 degrees. Attenuation and isolation as a function of rain 
rate were computed using both the- first-order multiple scat- 
tering and complete multiple scattering models. The results 
are plotted in Pig. 3 for 30 GHz. Note that there are only 
slight differences between .the values computed for the two 
models. At lower frequencies the differences were negligi- 
ble. From these calculations it can be concluded that high- 
er-order multiple scattering effects in rain begin to appear 
in the range of 30 GHz. However, results using the single 
scattering formulation of Section 5.2 would deviate notice- 


ATTENUATIO 
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30 GHz 


L=5 km 
uniform rain 
e=45° 

45° linear polarization 

< 0 >= 0 ° 


Og = 12° 


60% oblate drops 



— Multiple scattering 
• First-order multiple scattering 


140 


RAIN RATE (mm/ hr) 


Figure 3a. Comparison of multiple scattering and first-order multiple scattering rain attenuation 
calculations at 30 GHz. ' 
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ISOLATION (dB) 


Multiple scattering 
First-order multiple scattering 


30 GHz 


L-5 km 
uniform rain 
€-45° 

45° linear polarization 

<0>=O° 


00 = 12 ° 


60% oblate drops 


RAIN RATE (mm /hr) 


Figure 3b. 


Comparison of imltiole scattering and first-order multiple scattering rain-induced 
isolation calculations at 30 GHz. 




'-Vi. , r. S , 




. To further justify the claim that higher-order multiple 
scattering effects of rain are negligible. on communication 
links operating below 30 GHz, the calculated attenuation and 
isolation are plotted in Fig. 4 for parameters matching 
those of the VPI&SU COMSTAR D2 satellite, earth terminal. 

To more realistically represent rain, the piecewise uniform 
rain distribution of (34) was used. Mote the nearly identi- 
cal results with and without higher-order multiple scatter- 
ing included. 


5.5 Conclusions 7 : . ' 

In this paper formulations were presented for computing 
the effects of electromagnetic scattering from tenuous par- 
ticulate media. The formulations were for single, first- 
order multiple, and complete multiple scattering. Each 
included the distributions of particle si 2 es, shapes, and 
orientation angles. 

The multiple scattering formulation as presented here 
offers several features for calculation of millimeter wave 


propagation through precipitation media. These features are 
as follows, (a) The formulation is vector in nature and 


easily accomodates arbitrary polarization states for the 
input wave and receiving antenna, (b) Complete medium depo- 
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20 GHz 


L=5657 meters 

Piecewise uniform rain rate model 
6 s 44 ° 

Linear polarization 35° from vertical 

< 0 >» 0 ° 

cr =12° 

U 

60% oblate drops 


Multiple scattering 
First-order multiple scattering 


RAIN RATE (mm/ hr) 


Comparison of multiple scattering and first-order multiple scattering rain attenuation 
calculations for the Comstar satellite at 20 GHz. 







ISOLATION (dB) 



RAIN RATE (mm/hr) 

Figure 4b. Gonparison of multiple scattering and first-order multiple scattering rain-induced 
isolation calculations for the Comstar satellite at 20 GHz. 



ORIGINAL PAGE IS 
OF POOR QUALITY 



larization effects can be calculated (i.e., attenuation, 
isolation, and phase shift), (c) Scattering particle dis- 
tributions of particle size, shape, and orientation angle 
are directly included into the model, (d) Varying medium 
density along the propagation path (such as rain rate) can 
be accomodated, (e) Ice as well as rain hydrometeors can be 
included, (f ) The multiple scattering formulation is numer- 
ically efficient. For typical communication link calcula- 
tions the computer solution executes at least 30 times 
faster than the first-order multiple scattering version. 

(g) The multiple scattering model will give accurate results 
for systems operating above 30 GHz. 

The multiple scattering model was used to calculate com- 
munication link performance for rain along earth-space 
paths. It was shown that higher-order multiple scattering 
effects only begin to become important in the 30 GHz fre- 
quency range. 
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Chapter 6 

APPENDIX: COMPUTER PROGRAMS 


In this section three computer programs are presented. 
They correspond to the following: the simple , attenuation 

model of Chapter 3, the single rain drop scattering coeffi- 
cient derivation of Chapter 4, and the multiple scattering 
model for rain of Chapter 5. 
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C** ******************************************************** ************* 

c* * 


C* SIMPLE ATTENUATION MODEL * 

C* * 


c*** ************************************************************** ****** 


c* 



* 

c* 

THIS PROGRAM CALCULATES THE SLANT-PATH ATTENUATION DUE TO 

* 

c* 

A POINT RAIN RATE AT AN EARTH STATION. 


* 

c* 



* 

c* 

THIS PROGRAM WAS DEVELOPED BY THE SATELLITE COMMUNICATIONS 

* 

c* 

GROUP AT VIRGINIA POLYTECHNIC INSTITUTE AND STATE UNIVERSITY 

* 

c* 

AND WAS SUPPORTED BY NASA UNDER JPL CONTRACT 

NO. 955954 

* 

c* 



* 

c* 

ASSUMPTIONS: SPECIFIC ATTENUATION FROM OLSEN, 

ROGERS, 

* 

c* 

AND HODGE (1978). 


* 

c* 

POINT RAIN RATE DISTRIBUTION FROM CCIR 

* 

c* 

DOC. 5/5049-E (DRAFT RPT. 563-1, 

1981) 

* 

c* 



* 

c* 

INPUT PARAMETERS: 


* 

c* 



* 

c* 

FREQ = FREQUENCY (GHZ.) 


* 

c* 

ELEV = ELEVATION ANGLE OF SLANT PATH 


* 

c* 

LAT = LATITUDE OF EARTH STATION 


* 

C* 

HO = ALTITUDE OF EARTH STATION (METERS) 


* 

c* 

REGION - RAIN CLIMATE REGION 


* 

C* 

1=A 5=E 9=J 13=N 


* 

c* 

2=B 6=F 10=K 14=P 


* 

c* 

3=C 7=G 11=L ‘ 


* 

c* 

4=D 8=H 12-M 


* 

c* 

LABEL = OPTIONAL ALPHANUMERIC STRING (22 

CHAR. MAX.) 

* 

c* 

* 


c* *********************************************************** *********** 


c 

DIMENSION T(7),RAIN(7,14),LABEL(6) 

REAL LAT, L 
INTEGER REGION 

C* INITIALISE RAINRATE EXCEEDANCE DATA (FROM CCIR DOC. 5/5049-E TABLE I) 
DATA T/.OOl, .003, .01, .03,0.1,0.3,1.0/ 

DATA RAIN/22.,14.,8.,5.,2.,1.,0., 

$ 32. ,21. ,12. ,6. ,3. ,2. , 1. , 

$ 42 . , 26 . , 15 . , 9 . , 5 . , 3 . , 0 . , 

$ 42 . , 29 . , 19 . , 13 . , 8 - , 5 . , 3 . , 

$ 70. , 41. , 22. ,12. ,6., 3. ,1., 

$ 78., 54., 28., 15., 8., 4., 2., 

$ 65., 45., 30. ,20., 12., 7., 0., 

$ 83. , 55. , 32. ,18. ,10., 4. ,0.', 

$ 55., 45. ,35. ,28. ,20. ,13. ,0. , 

$ 100. ,70. ,42. ,23. ,12. ,6. ,2., 

$ 150., 105. , 60., 33., 15., 7., 0., 

$ 120., 95. ,63. ,40. ,22. ,11. ,4. , 

$ 180., 140., 95., 65., 35., 15., 5., 

$ 250 . , 200 . , 145 . , 105 . , 65 . , 34 . , 12 . / 

C* READ THE INPUT DATA VALUES 

READ (5,5) EREQ,ELEV, LAT, HO, REGION, LABEL 
5 FORMAT (4F10 . 3 , IX, 12, 7X, 5A4, A2 ) 
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C* CALCULATE SPECIFIC ATTENUATION COEFFICIENTS 

IF(FREQ.GE. 1. -AND. FREQ. LE. 1000. ) GO TO 20 
WRITE(6,10) FREQ 

10 FORMAT ( 1H1 , ' FREQUENCY=* , F10 . 4, * IS OUT OF PERMITTED RANGE*, 

$ ' - - PROGRAM STOP*) 

GO TO 999 

20 Al-6. 39E-9*FREQ**2 . 03 

IF ( FREQ . GE , 2 . 9 ) Al“4 . 21E-5 *FREQ**2 . 42 
IF( EREQ . GE . 54 . 0 ) A1=0 . 0409*FREQ**Q . 699 
IF(FREQ.GE. ISO. ) Al=3.3S*FREQ**(-0. 151) 

B1=0. 851*FREQ**0. 158 

IF(FREQ.GE.8.5) B1=1.41*FREQ**( -0.0779) 

IF( FREQ. GE. 25.0) Bl=2. 63 *FREQ**( -0.272) 

IF ( FREQ . GE . 164 . 0 ) Bl-0 . 616*FREQ**0 . 0126 
C* INITIALISE RAIN DISTRIBUTION PARAMETERS 
GAMMA=l./22.0 

C CALCULATE AVERAGE 2ER0 -DEGREE ISOTHERM HEIGHT 
HI-4.8 

IF(LAT.GT.30. ) HI=7 . 8-0 . 1*LAT 
C SUBTRACT EARTH STATION ELEVATION 
HI=HI- (HO/IOOO. O) 

ELEVR=ELEV*3 . 14159265/180. 

C* PRINT OUTPUT HEADER 
WRITE(6,25) LABEL 

25 FORMAT ( 1H1,///, 11X, ' SAM - - SIMPLE ATTENUATION M0DEL'// / 5X,5A4,A2) 
WRITE(6, 30) FREQ, ELEV, LAT, HO, REGION 
30 FORMAT (' O', 4X, 'FREQUENCY: ',F5.2, f GHS f // 5X, 'ELEVATION: ',F5.2, 

$’ DEGREES ’/, 5X, ' LATITUDE: ' ,F5.2,' DEGREES’,/, 

$5X, ' ALTITUDE: ',F5.1,' METERS' /5X, * RAIN CLIMATE REGION: ',12) 
WRITE(6,35) A1,B1 

35 FORMAT( 'O’ ,4X, 'SPECIFIC ATTENUATION COEFFICIENTS: A2=',F8.5,/, 
$40X, T B1=' , F8 . 5 ) 

WRITE( 6, 40) 

40 FORMAT ( 1H0 , 7X, ' PERCENT ' , /, 9X, ' TIME RAINRATE ATTENUATION 1 / ) 

C* CALCULATE ATTENUATION 
DO 70 J-1,7 
N=S-J 

R=RAIN(N, REGION) 

HE=HI 

IF(R. GT. 10. ) HE=HI +ALOG10 ( R/10 . ) 

L=HE/S IN ( ELEVR ) 

ALPHA=A1*R**B1 
IF(R. GT. 10. ) GO TO 45 
A=ALPHA*L 
GO TO 50 

45 ARG=GAMMA*B1*AL0G(R/10. )*COS(ELEVR) 

A=ALPHA* ( 1 . -EXP ( -ARG*L) )/ARG 
50 WRITE (6, 55) T(N),R,A 

55 FORMAT (9X,F5. 3, 5X,F5.1, 7X,F6.2) 

70 CONTINUE 
999 STOP 
END 
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6.2 SINGLE RAINDROP SCATTERING PROGRAM 


The COEFF program was written based on the algorithm 
developed in Section 4.2. It computes the scattering coef- 
ficients of a raindrop. In this section we will describe 
the main subroutines of the program. 


i r 


Subroutine GENLGP. 

This subroutine calculates the associated Len- 
gendre polynomials. It uses the recurrence formu- 


5^(x) = 1.0 P°(x) = x 


(2ti-l) I? , (x) x - (irim-l) i? 
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published in literature. (1) First we compute the P 
and Q polynomials# then we compute the spherical 

7 

Bessel functions and h£(x). P and Q are given 

PY 

[1/2 n] u _ 2 k 

P(n + 1/2, x) = £ (-l)^(n + 1/2, fc) (2x) ^ (5) 

fc =0 


[1/2 (n-1)] -2k-l 

Q(n + 1/2, x) = £ <-l) K (a + 1/2, 2k + 1) (2x) ^ L (6) 

where 

(n + 1/2, k) =• —fc i) 

Then j n (x), y n (x) are given by 

j n (x) = x“^[P(n + 1/2, x) sin(x 5p) + Q(n + 1/2, x) cos (x~?f )] (8) 



( 9 ) 

y n (x) = (-l) 1 **^ x”^[P(n + 1/2, x) cos(x + 5JL) - Q(n + 1/2, x) sin(x + ?p)] 


(l)Abramowitz, M.# Stegun, X.A., Handbook of Mathematical 
Functions, Doven Publications, New York, 1965. 




Then 


h n ^ = * j y n (x ^ 


3. Subroutine GENGEN. 

This subroutine calculates Gegenbauer polynomi- 
als using the recurrence formula 
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( 10 ) 


<$(x) = 1.0 • .. 


( 11 ) 


C^(x) = 


2. * x 


( 12 ) 


= < 1 < X > * 4-l (x) - 


(13) 


4. Subroutine VM&TR. 

■+ ■+ 

This subroutine calculates VCfe^# ^ 2 ) as 


sinlfc - So| - cos|l£t - So | 

V<kp - 4irabc y — — ~ T 2 ' (14) 

I*L " 


Note when = £3 we divide by 0 so we must 

■> 

% take the limit &2 k l giving 


"3* / 

UO-c^, kj) = ~ tt abc y 


(15) 
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5 . Subroutine ZMATR. 

This subroutine calculates the integrand of 
(27). 

6. Subroutine KNTGR. 

This routine performs a numerical integration 
over x in (27) using a Clensaw-Curtis quadrature. 
Then it calculates matrix K(ki, fc 2 ) 

K = s r I V<^, %) - l<jt, %) ( 16 ) 

The seven point Clensahw-Curtis quadrature is 
given as; 

f 1 7 

/ f(x)ds=S c f(x ) (17) 

J -1 n=l n “ 

where: 

c L - 1/35 = c 7 c 2 - c 6 - 16/63 c 3 - c 5 = 16/35 c 4 -164/135( 18) 


- -1 , X2 “ -1/3/2 , Xg ” -0.5 , x 4 = 0.0 , Xg = 0.5 , (18) 

Xg = /J/2 , Xj - l 

Xn the main program we develop the matrix equations to 
solve for [C(1ci)3 and then calculate for the scattering 
tensor £(■& , ) . 




Inputs to the program ares the equi volumetric radius of 
the raindrop, the length of the major axis, frequency of 
incident wave , relative dielectric constant of raindrop, and 
the directions of the incident and scattering waves. The 
output is the scattering matrix of the rain drop. The list- 
ing of the program follows. 



TSll >0o2 
IS" 000 J 
l STi 0 0 u h 

isw on os 
iso uuoo 

i Sf! 0007 
ISO 0(< Jd 
ISO OC'tip 
IS** 0 tj 1 0 


CSJOn 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

C I ms PROGRAM CALCULATES THE SCAT t RING it OS OR OF A D] ELECTRIC 
C SPHEROID. THE VARIABLES ARPF. AS FOLLOWS: 

c 

C THE I , PHI 1 : THE DIRECTIONS OF THE INCIDENT WAVE 

C THKS,PrilS : 1HF D I KECTtCMS OF Thfc SCATERED WAVE 

C AIJAhSTME equivolhmktric RADIUS OF THE SHFPOID 
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C T»i r . PhuGPA 1 * * ’ WAS WRITTEN R i A? A51 ASICS TSOLAKIS. 

C THIS PROJECT ..AS SUPERVISE? RY TR. V.. L. STUTZMAN 
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C 

c 


: HI. f C I 1 KEAL*8(A-h,E-H,U-?) .COKPI LX*16(C) 
C ' • * 4 P L F X OK K (60, 60 ) , CUN4 T ( 60 , 3 ) 

,J ~ A L * 1 o FACT ( '20) ,FACT1 

;3-:al nwA(60) 


I) f HE4& J ON C*MATR( 3, 3) ,CAJ(3.3) 
C#CF(3,i) # AFACT(20) ,CSUMl3#3) , A J ( 3 , 3) 
C'HitfG.J/ACi/ AFACT .FACT 
C IMON / if- 1 1 G H / C’i ( 7 ) , W C 7 ) 

C VtMns/PrtH/ A # A C , C S A * • , P I # C F P S 1 1 
P I = i . 14)502*353560 7 7 
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/ 1 = 1 

4 = 9 

he:;; ite varamlfs 

'• ap (5.90) AhAH.A.FO.CI R 
0) F P“tl (2FH.6.5X.r0.2,7(5X,Eb.3)) 

• ITc(n f H9 ) ABAP, A # FO,CNO 
A«»(5,lo] ) PmT J , THE 1, PHIS, Ihfb 
»<; A U4f I 0.5 ) 

(Tf l ", 10 | ) PHI 1 , THF I ,PHlS,ll Ff 
■ L* »C = 3 . H-0 l 


ORIGINAL PAGE t& 
OF POOR QUALITY 


Cm 3 


*• ^ ^ ! ! * ' C? ! ! l Cm3 C3 CmTI CT!^I ITTD Cm3 




00 

1 


ISN 00?'?. 
TSN 0023 
ISli 0024 
ISO 0025 
ISO 0026 
ISO 0027 
ISO 0028 
ISN 0020 
ISO 0030 
ISO 0031 
ISO 0 u 3 2 
ISN 0033 
ISN Ou 34 
ISN 0035 
ISN 0036 
ISO 0037 
TSN 0036 
ISN 0036 
IS i 0040 
ISN 00 Hi 
ISN 0042 
ISfi 00 h 3 
ISO 0041 
IS!! 00 4 6 
ISN 00 4 0 
ISO 0047 
ISO 0048 
IS 4 00*9 
ISN oob) 

150 00S1 
IS.J 0062 
IS!, 0063 

15 1 Gy 6 4 
IS* 1 0056 
ISN 0 0 6 o 


ISN 0057 
131 0058 
J S .i 0 059 


IS 1 1 OOoO 
IS'j 0061 
1 S ( 0 ( i r> 2 

1 .5 i tj*>o ! 

1 ., i 0 1 1 <> ; 

I S • ,» Mr S 
ISN 0 f i o • ' 
1 S i l 1 9 1 
l s •: o o f- >* 

T 5 i *»«>»» y 
IS 1 i’ WO 


V 

C 


n 


f r=(ABAR**3)/(A>M2) 

/■C« AC/100. 

A = A/ 1 00 . 

A ‘■'0 = 2 .«P I *H)/tfFLOC 
PHII=PHII*PI/180. 

THE 1=1 H£I *PI / 1 80 . 

OHIS=PHlS*PI /180. 

1 H£S=lHES*PI/lbO. 

O:PSlL=Cf*0**2 
O l *Tfc:i=3.0D + 00 

h. SQ=[>SQP I (RTMK.I ) 

.V E = 1 .01)400 

C.r, NL=f 1 .Q.O.O) 

YY( Us-ONE 
i . (3)=-G.vE/2. 

vv( 2 )=-o!Jt*psg/ 2 . 

i V( 4 )=QNF.*0.0 

i . (5)=-YYUi 
yv(b)=-V *12) 

Vt(7)s(JLE 

r ’( l)=C0NE/35. 

c 7 ) =C^- ( 1 ) 

C . (2)=CUiJE*l b/6 3. 
r-.m=C*(2) 

C. ( 3)=CCNF*16,/35. 
r ’ ( 5 ) =C ». ( 3 J 
C ( n=C0wF.*lo*./J15 t 
C',A 4=< AF U**2)*(CEPSjL-l . > 

C A K l = A h 0 *C N Q 

f AK 2 = AKG*Cwn 

f? AKOsCONh* AKO 

"i=2*v*2 

A” 5 1=1*02 

»• cm )=r am u-i ) 

5 i- ACT I I )=Ar ACTl ( I-l) 

CALCULATE IHK JAI TFNSOR IN THE PT SECTION GF THE INCIDENT 

.v A 7 F 

CALL A J v A ( PH 1 1 « 1 HE I # A J ) 

T TGU W=-3 
ir )":.2 = -3 


KMAbLlSh Ittt: -* A T H I X EQUATIONS lhAT THE 1 *.<0 COUFLFP INTEGRAL 
CP- 5 P L E D EQUATIONS HAVE BFFiJ PEDUCFD 
T N r DT-lENSlON UF THE DKK MATi TY IS NX*NPHM9*3 

• I J 1 = 1 t NPh 

* -1=PTM 1-1 . )/(NPI,-l. ) 

•y h 1 J = 1 /NX 

< I =-l . + 1-1-1 . )*2./(NX-l . 1 
l .'JU;.l =J CHUM 4 3 

* 20 f\ = l,HPH 

* '2 = 1-14 lr -1 . )/(LPH-l . ) 

• u l.= 1 , rlX 

» •> = •!. 4(L-1. >*2./(«*X-l.) 

I '.i.i, 2= I Ci li; i ? 4 3 

WGPlf.XI ,CAK1 ,CAK? ,PH1 # l P2,X1 , X ? #C AhO , N ,CKMATR ) 


ORIGINAL PAGE IS 
OF POOR QUALITY 




r 


v£> 

I 


ISM 

0071 


L.i 3 0 1 r = 1 , 3 

IS'J 

0072 


DP 30 1 J = 1 ,3 

isn 

00 7 3 

30 

C 30 

OKKUCDDM + IK , TC0UM2+1J)=CKMATM IK ,IJ) 
DKM 1 CUUN ? + 1 i*v , j COON 1 + 1 J ) sCKPATR ( 1 h , IJ) 

ISH 

0074 

20 

cnNTinut: 

ISN 

0075 


■vR r TE ( 6 , 600 ) 

ISN 

0076 

6 00 

F ' i K M A T ( 5 X , 5 H K N T G H ) 

ISM 

0077 


lCUUN2=-3 

I si: 

0078 


T'u =dakcps(xi ) 

TSiJ 

007 9 


CAM. U A I K ( C A K i , P H 1 , ’1 H 1 ,CAKO,PHJ 1 ,'IHEI 

ISN 

CO 60 


CM. 1. CMIM.T ( A J , 3 , 3,CM,CAJ) 

T S -1 

008 1 


of’ 5 0 IK 1 = 1, 3 

1 Si! 

0 0 b 2 


.H. 50 l J 1 = 1 , 3 

OHMITf 1COU.4 1 + IKI , MI)=CAJ(1KI ,1JI ) 

I S.J 

00 b 3 

50 

I so 

00 84 

10 

CiJN ril'Ut 

ISH 

008b 

c 

lMs;JX*NPH*3 


ISH 

iso 

I S -4 

is g 
ISM 
l Si; 
ISM 
I 6 i 
is*: 
1 Sti 
1 SM 
ISn 
1ST 
ISM 
I S'! 
i S‘i 
I s 'I 
I Si 
J SM 
1 S -i 
IS, I 
) Si I 

IS i 
l S : > 

1 S. i 
I S * 
I S I 
TS i 
1 Sf. 

( s • 


0086 
0087 
00 8 8 
0 0 H 9 
0090 

009 1 
0 0 9 2 
3o9 I 
00 9 1 
0095 
0098 
00 9 1 
0 0 9 8 
0 0 9 9 
0 l G 0 

010 1 
010 5 

0 1 0 j 

u 1 04 

OIOS 
0 1 <• c 
0 1 1- / 
0 1 o < 

0 1 K * 
0 1 1 i: 

Mill 

o 1 w 

0 1 i I 

■ i 1 j -• 

O 1 1 -> 


C SOf-vE FH8 I HE C ( K ) MATRIX AND THE I ! CAICUI.ATL THE CF(3,3) 
C SCAfEKUG ILUSOR in THE DIRECTION OF THE SCATEHED WAVF 

'.'RITE! 6 . 700) 

700 f n RMAT(8X , 6 H L EQ T 1 C ) 

CaLI. L£tiTlC(DKK,r|fc f NN,DUnl IT ,3 , 3, 0,08. A, IFF) 

' 6 I TE ( b , 126) I FR 
128 c ' r 'P M A7 (5X t *THK TER 1S # ,I5) 


c*5 


HO 


7 0 


ion 

1 o 2 


•? I j t. ( b , 7 0 0 ) 

D" 65 1=1,3 
r.- 6b J=1 , 3 
C 3 Li M ( J , J ) = 1 0 • 0 , 0 • 0 ) 
i r r jM'i=-i 
i : 7 o 1 = 1 , Mpi-j 

i =Ri * ( J - i . ) /(N.Pn-1 . ) 

LG 7 0 J = 1 , N X 

VI =-l .+ (0-1 . )*2./(*iX-l. ) 

HO K=l,3 
V) 8 0 L= 1 , 3 

'V* J ( K , L ) =Dtl.'JI T ( TCOUfi + K , L) 

1 RE1 sDAhCL»S( / 1 ) 

C - 1 * L Iv-iATR ICAKfi, [■ H J.S , THtS ,C AK 1 , FH1 , 1HE1 ,CU) 
C.^IiLf C u ,k Uj L 1 (cAJ,3,3,Cl.) 
call CSU 1 -' (C.smh, :i ,;.|,CAJ1 

C-> I 1. AwDiA 1 Ril I S , 1 .if.S , -*J ) 

i i. iOl.UA.I, i, 3,CA ), 3,CK) 

ICO 1 = 1,3 
):' 100 J=l,3 

■ *M IT ( o, 1 o 2 ) CF( i ,1) 
c ' i3 ) A V ( 5 X f 3 ( t 1 2 . 1 , 3 X , c! 1 2 . 4 , 2 X ) ) 

■» t V.’P 

•. ' t . 


ORIGINAL PAGS IS 
OF POOR QUALITY 


CZH' =3 EZ 3 1=3 1 = HZ 3 r=I £= 1 = =1 m 




ISN 0002 


|S.U 0003 
ISM 0004 
ISM 0005 

ISN 0006 


ISN 
1 5 M 
ISM 
i*i!i 
TSH 
ISU 
ISM 
ISM 
ISN 
ISN 
ISM 
ISM 
ISN 
ISN 
ISM 
ISN 
ISM 
ISM 


IS:* 
IS!' 
IS "4 
ISN 
IS i 
IS’J 
ISM 

IS -J 
IS* 


IS 

IS 

IS 


00 0 7 
0008 

0009 

0010 
0 011 
0012 

0013 

0014 

0015 

0016 
0018 
0019 
0021 
0022 

0024 

0025 

0026 


00 2H 

0029 

0030 

0031 

0032 
00 33 
0031 

00 36 
003o 


IS* 0037 


00 3 w 
0 0 39 
U 0 ] 


C 

C 

C 


C 

C 

c 

c 


SUBROUTINE ZMATRXCPH1,X1,PH2,X2,CAKI , C AK2 ,CAKO, X , N ,CZ ) 

THIS SUBROUTINE CALCULATES THE 21 MATRIX 

IMPLICIT REAL*B l A-B,D-H, 0-2) , COMPLEX*! 6(C) 

HEAL*!* FACT f 20) 

DIMENSION CZC3,3) ,CAK1J(10) ,CAH 2J( 10) ,CYK0Y ( 10) ,CYKGJ( 10) 
C #CSUM ( 3 * 3 ) 

DIMENSION SI)M(3,3),PCHC10,10),PCHH10,1Q),PCH2(10,10) 

C ,GEN C 10) ,0(3 f 3) , AFACTC20) 

COMMON/ ACT/ AFACT.FACT 

COMMON /PAR/ A j AC,CGAM, PI ,CEPSIL 

t\IAlsfC.O r i s O) 

S'i-A**/ 

/CSAsAC+*2-A**2 
nn 5 1 * 51,3 
no 5 1 J=1 ,3 
I CZ(IK,IJ)=(G. 0,0.0) 

ACH=AC*X/P$0RT(SA+ACSA*eX**2)) 

1FCDASS t ACH ) «EO* 1 « 0) RETURN 
aCHI=aC+X1/DS0RT(SA+ACSA*XX1**2)) 

I. F ( D APS ( AC H 1 ) ® E 0 • 1.0) RETURN 
ACH2~AC*X2/nSQPT(SA+ACSA*(X2**2)*) 


C 

C 

c 

c 

c 

c 

c 


TFCDABSCACH23.E0.1.O) RETURN 
Cf»*Kl=CAKl*DSQRT(SA*C 1 •-Xl 3 * t *2) + (AC*Xl )**2) 

2 ** 2 )+ 

X)**2) 


CAKK2=CAK2*DSQRT(SA.*(1 --X2**2) + (AC*X2)**2) 
YsDSQRl (SA*(1.-X**2) + (AC* 


AND THF. SHERI CAL 


C YKU-C AKu* ¥ 

CM.CULAit* IMF LERGENORE FUNCTION 
SSFL. FO lCTItJUS OF ORDER N+l 

t | , - $ ^ \ 

CALL GEMGP ( PCh , ACH, NN ) 

CALL GEFLGPf PCril, ACH1 ,NN) 

CALL JEN LGP C PC H 2 , AC H 2 , N k ) 

CALL SHPSESCCAKK1 ,CAK1 J ffCYKOV ) 
r 'LL SHPBF:sCCA^'t<2 f M^:CAK2J,CYKiY) 

CALL SRP6ESCCYh0,MM,CYK0J tf CYT0Y) 

!t PlT£f 6,300) UN 

300 FORMAT (5 XV’TriE ORDER IS ',13) 

\r U2=PH 1-PH2 

CALCULATE i HE 3t£ iGt 1 * R ATFEH POLYLPMaLS 

r „LL GENGEJJ c Or ** , DCOSC PH 1 2 ) , N + 1 ) 

TkVCATF THE INFINITE SUMP ATI Oh TC GRRFR N AND DO 
«L' THE Cii .PHTAUOi-iS T>1 FIMp THE 7 1 MATRIX 

I s V, * l 

1 • 1 \» I S 1 , M l 

' f It) 0 “ 1 , J I 


ORIGINAL PAGE 
OF POOR QUALITY 


t 


ISn 0041 
I s :i 0042 
T SO 004 3 
ISO 0045 
TSIJ 00 4b 
ISN 0047 
ISO 0048 
ISi'J 00 49 
ISM 0050 
ISO 0051 
ISM 0052 

ISfJ 0053 

I Sli 0054 
1ST 005b 
ISO 0053 
1 S i u 0 5 9 


0060 
0062 
0 06 3 
0 06 4 
0065 
00 66 
0067 
OOOH 
0 0 70 
00 7 1 
0 0 7 2 
00/3 
007 4 
00 75 
0 0 7 c 
00 7 7 
00 7 8 
0 0 79 
CO 60 
OOfe 1 


I :vi 
1 S j 

J SO 
1 o 1 
I S': 

1 O. 

1 i.. 

IS' 

ISM 

J S . 
1 6 
I s* 

lzd • 


0 08? 
0063 
0 0 a 4 

0 n ft 5 

* 1 1 .) 4 n 

< < 0 h / 

0 1 ) ►* 4 

0 n 9 
0 09 »j 
l* 09 1 
O'# 9 2 
P i - 9 1 

GD* 


i, i=I-i 
1 1=J-1 

TF(I(IN+lN)/2)*2..4E.IWt!«) GO 10 10 
no 20 16=1,3 
03 20 I J= 1 , 3 
) SUM (IK,IJ)=0.0 

on 30 K= 1 , I 
0*1 30 L=1,J 
XS=K-1 
I.T = L- l 

ACQU1 = ( K S + 1 )*AFACT ( 14IN-KS)* (LI +1 ) *AFACT ( 1 + IM-LT ) / 

C( t FACT l iU + KS+3 J )»AFACr(IMtLTf3) ) 

ACON2=ACONl *PCH(IN+1 ,KS+1 )*PCH(IM4 1 , LT+ 1 ) 

C * ! J C H 1 (Jlu4l,KS4l ) *PCH2 ( IM+ 1 , HT 4 1 ) 
l F ( K S . Gfe*. • LT ) PHB=PH1 
T *• (KS.L1 ,LT ) PHd = PH2 
I •iI.J=hT?. ! 0(KS,l/f) 

G=2. H'MGEN C IM'H 1 1 

EVM.l'A’It THE II, 12,..., 16 INTEGRALS DEPEHDIG ON THE SUM KS , LT 
tt£f!G EVEN OR ODD 

IFC ( (KStLl J/2) ♦2.NE. (KS+KT)) GO TO 100 
All =G 
ft I 2 = 0.0 
A I 3 = 0.0 

1 4=G» (£>S1MPHB1*»2) 

ft l 5 =G* US i 0 ( PHh ) ♦ OCUS ( PUB ) 

A16 = G»(0C(JS(PHB)»*2) 

TF(KS.NE.LT) GO TO 200 

hI4 = A1‘* + PI *DC0S(KS*PH1242.»FH1 ) 

A15 = A I5-FI *081. HKS^Ph 124 2 . ♦Pin 
l I 6 = Aio-FI*UCLlSCKS + PHI24 2.*r hi) 

-«■ ro 20u 

1 I 1=0.0 

ft ! 4=0.0 

\ J 5 = 0 . 0 
A 1 6 = U . 0 

*M? = G4£»SIL(PHrt) 

4 ( 3 = iS*OCOS(PHH) 

I S V = 1 - X * ♦ 2 

.>'.-RY=DS0Rl (Si J 

CMCUHAlfc 1 HE O(S.T) MATRIX 

»l 1 , 1)=AI 1 - S Y * A 1 6 
•f 2 , 1 )=-SV ♦Alb 
.?( 3, t )=-X»SQRY*AI3 
-( l , ?)=-SY*M5 

' ( 2 , 2 j = A l J - S i’ ♦ A i j 

• f 3 Syr Y* A I? 

f 1 , 31=-X*.SfVj * A 1 3 

•• < 2 , 3)=->*S(.LY*AI2 

• f 3, 3)=Sl*All 

.‘ALL ft.UL l C ( Q , 3 , 3 , AGU'12 ) * 

' * U j. SUvCSHO, 3,1,0) 

~ iVl I M E 

j i tZM C3DMEED ED ZD ED LTD 


o o 
-n 2 
-o © 

os 

O ’iy 

73 r* 

o -o 
c £ 
> © 
r * rs 



r— i 


rssi rr-i 






1 SN 009b 
1SN 009o 

1 5N 0U97 
ISM 00 9 3 
ISM 0099 
ISO 0100 
ISM 0101 
ISIJ 010? 
I SiM 0 1 ft 3 


T "• I h = M I [JO ( lMalN) + ) 

C'0|ils(2. + Iflf 3. )*12.*IM+3)*(CAK1 J ( I N+ 1 ) /C AKK 1 )♦ 

C(CAK2J( I MM )/CAKK2)*CYKOJ(lMAX)*(CYKOJ(IMlN)+CJAI*CYKOY(lMIN)) 
CALL C f *iULT ( SUM . 3 • 3 , CON 1 , CSUM ) 

CALL MCSUMCCZ, 5,3, CSUM) 

10 continue 

Cnf;Y=CJAI/(Y**2) 

C LL CChULT(CZ, 3, 3 f CONY) 

.*FTlJRrC 
. ■ L> 










ISM 

000 2 

I S h 

0 o 0 3 

I so 

000 4 

1 SM 

0 0 OS 

T Si. 

ooos 

L Si 

0 0 o 7 

1 ‘v. 

OOUb 

I vl 

0 004 

ISO 

00 10 

1 S-‘i 

001 1 

IS J 

00 13 

1 SI! 

00 14 

IS I 

■T 0 1 5 

IS.J 

•JO is 

T Si 

0017 

I SI. 

0019 

LSI 

J0 2 l 

IS.. 

00 2 2 

r s 

00 2 3 

I Si* 

0 o 7 -* 

T Sr. 

0 0 2 S 

J 3*. 

00 2r> 


C 

c 

c 


b r ( I # ))—0«0 

irCASiL.EC.O.Gj RETUPN 
P( 1 , 1 ) = 1 .5/ AS IN 
P ( 2 , 1) = X / A S 1 M 

1 »«'< 10 i = 2 , t'j 

» / • 10 J = 1 , I 

1 K( I .EC* 2. ANP.J.EG. 1 ) GU TO 10 
IFCl.NK.J) GO TO 20 

M 1 . n=AFACT(l+2»(I-l ))*AS1.N* ♦CI-2)/( (2**(I-1 ) )*AFACT(I ) ) 

v*-' ru i o 

20 *'(l ,.l) = 1(2. ♦ 1 -3 ) +A* P( 1-1 ,J)-( 1+0-3 )*P(1-2,J)) 

C/i l-.J) 

10 C t*lL 1 I mUE. 

•M r; | Up\ 

~:.n 


subroutine geulgp(p,x,n) 

THTS SUhBUUTINE CALCULATES THE LENGF.MDPE FUNCTIONS 

I. PL1C1T KEAL*rt ( A-H , 0-Z ) 

REAL ♦ 1 6 FACT ( 20 ) 

.•MF USiriM P(M,U) . AFAOC20) 

CiMMOn/ACT/ AFACi $ F ACT 
ASlNsOSOfll (1 • -X**2) 

U n 5 1 = l # N 
*'•: s j = i,n 


o ° 

^ 50 


TJ 

O 

O 

X> 

o 

c 

> 

r~ 


2 

7t 

r- 

S 

O 

FI 


5 CO 




















JSN 009b 

isn 009o 

13N 0097 
ISM 0098 
JSN 0099 
ISM 0100 
ISM 0101 
ISU 0102 
ISN 0103 


10 


T*-'.lNs.yiIHO(IM c lN) + I 

Cniil = (2,+IfJ + 3.)*(2.«IM+3)*CCAK1 J(IN+1 )/CAKKl }* 

C f CAK2J ( I 1 ) /CAKK2 ) *CYKGJ (IMAX ) * (CYK0J( IMIN ) +CJAI *C YKOY ( IMIN) ) 
CALL C''iULT(SUM.3.3,CaNl,CSUM) 

TALL MCSUM tCZ , 3 , 3 ,CSUM) 

CONTINUE 

CnNy=CJAI/(Y**2) 

CALL CCfcULTCCZ , 3 , 3 p CONY ) 

.iFTURU 

**■■0 


ISM 0002 


I Sly 0003 

150 000 4 
lS f i 000b 

is*. 000b 
I S i 0007 

LS't 0 0 0 b 
I 3U 0 0 09 

iso 0010 

IS '"i 001 1 

151 0013 
ISU 0014 
ISI 0015 

15.1 0016 
IS.* 0017 
JSfi 0019 
13l 002 l 

13.1 0022 

rs.i 00 2 3 
I <S 0 u 2 •* 

isr* 002s 
JS'f C*G2n 


c 

c 

c 


SUBROUTINE GENLGPfP,X,N) 

THIS SUh RUNTIME CALCULATES THE LEMGENDRE FUNCTIONS 

f. PLicn real*hca-h,o-z) 

REAL *16 FACT (20) 
i.MFOSinfi p(U,!jJ , AFACIC20) 

CiiMnaw/ACT/ AFAC’i ,FACT 
ASlrJsnSCftl (1 .-X**2) 

L> f ' 5 I = l,N 
5 J = l,l. 
b rtl,j)=0,0 

LF(ASIL.KCi.O.O) RETURN 
R f 1 , ! ) = 1 .O/ASiN 
H C 2 # 1 ) - X / A S I ‘ J 

Mm 10 1 = 2, iu 
Li 10 J=1,I 

IF(J»EC.2.4Nn.J,fc;Q.1) GO TO 10 
IFll.NK.J) GO TO 20 

, l)sAFACT(l+2*(I-l ))*ASlN**(I-2)/(C2**£I-l))*AFACT(I)) 
fu 10 

20 ^Pf j ,J)=l (2.*1~3 )*a*P(!-1 ,J)-(I + U-3)*P(1-’2,J)) 

it> 'f<UlJ.vUE 
) l'uR'y 
-:\o 
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ISM 000/ 


1 1'Bm 


000 3 
OOG-l 
00 ob 
00 06 
000 7 

0 ')*') h 
0 0 0 0 
00 1 “ 
I»01 1 
0012 

001 3 
00 14 
0 0 l 5 
00 In 


00 17 
00 1H 
0 0 1 9 
00 2 n 
0021 
0022 


00 2 3 
0 0 24 

0 0 2 d 

o 0 2 n 
o02 7 
002 7 
t» u 3 1 
0 0 4 3 

00 . 3 4 
1 » 0 3 S 
0 0 3n 
'10 i I 
0 0 4 
00 4 9 

00 I*' 

0 0 4 l 
•' ! • * *♦ £ 

1. I * 4 

1 o ; -i 
1 < 0 4 •> 

"h , •-* 

• • 0 H / 

M • J In 

• • ’ * t ! 






t wm*** 


S 1 1 R K0UI i ii E S H PB FS ( 0 A RG ; N , 0 JSH , 0 Y SP ) 


THIS Slip ROUT | NE CALCULATES THE SPHERICAL BESSEL FUNCTIONS 
THE ilFTHOD USED IS THE ONE DISCRIUFD IN THE! ABRAMOV I TZ, STEGUN 
riA o>finoK OF MATHEMATICAL FUNCTIONS. QUADRATIC PRECISION IS 
MECESSAR* 


I 'PLICIT REAL *16lA-b,D-H,P-Z) , COMPLEX* 32 (C ) ,C0MPLEX*16(0) 
tKAL*R AFACT ( 20 ) 

DIMENSION OJSH(N) ,GYSH(N) ,FACT(20) 

C n -‘. r iOM/ ACT/ AFACT, FACT 
J/,HG=OARG 
vi 20 I = 1 ,M 

rrs( 0 . 0 , 0 . 0 ) 


.:o=co t o f o.o) 
? i = l / 2 ♦ I 


• 2=( 1-1 )/2M 
on 25 Js If 01 
4 ^ = J - 1 
-=J-1 

rHa((-l.)*FK)*FACI , (ltIf2*K)/(FACT(l+2*K) 
C * F AC1 (1 +1-2* K)7 
CP=CPtCPl/l (2.*CARG)*« (2*K) ) 

CONTINUE 

.)•' 30 J = l,.i2 

ansJ-1 

rs=J-l 


: 'l = Cf-l . M+Km , ACj , (I+2*K + 2)/(FACT(2*K + 2) 
C ♦FAC h J • 2 ♦ K ) ) 


w >=C0 + C 0 1 / ( ( 2 • +C AHG J ** ( 2 *K + 1 )) 

C».h I i MOL 
n ’ i = I / 4 . 

4=1/4. 

Tr'f AU.Ey.N 3 + 0.25) GO IU 60 

T *• C A ! l • I* v • *j 3 4-0 . 5 ) GO TO 7 o 

irl AH.EQ. >*3 + 0.75) GO TJ 80 

°>4P=lCP4COSll (CARG) +CU*CQCUS (CAKG ) )/CARG 
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6.3 MULTIPLE SCATTERING RAIN PROPAGATION PROGRAM 


The multiple scattering formulation for the rain problem 
as developed in Chapter 5 has been coded into Fortran. The 
block diagram of the Rain Multiple Scattering Program (RMP) 
is shown in Fig. 1. The correspnding analytical development 
and examples of computer results are discussed in Section 
5.3. 

The statement listing of the program follows. 



PLOT ISOLATION , 
ATTENUATION vs. R 


CALCULATE MEDIUM 
SCATTERING TENSOR 
fw 


SUBROUTINE 

COEFE 


FIND El GEN- VALUES , 
VECTORS OF K 


SUBROUTINE ' 
EIGVC 


SUBROUTINE 

MULT 


SUBROUTINE 

OUTANT 


CALCULATE 

ISOLATION, ATTENUATION 


SUBROUTINE 

MAP 




Figure 1. Block diagram of KMP program. 
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Input parameters of 
rain medium aud incident 


Calculate Ex, Ey 
components of incident 


SUBROUTINE 
COMP NT 


Output input parameters 


R - 10, 160, 10 
o 


R rain rata in tam/h 
o 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C THIS PROGRAM CALCULATES THE E VECTOR FIELD THAT COMES OUT 
C FROM A RAIN CELL. If ALSO CALCULATES THE POLARIZATION 

C DIRECTIONS OF THE INCIDENT WAVE THAT WILL PASS THROUGH 

C IHE RAIN CELL UNCHANGED. 

C El IS THE INCIDENT FIELD TO THE RAIN CELL 
C I F I V IS THE VECTOR FIELD COMING OUT OF THE RAIN CELL 

C L IS THE THICKNESS OF THE RAIN CELL IN METERS 

C LOKANG IS, THE ELEVATION ANGLE OF THE SATELLITE 
C FREQ IS THE FREQUENCY OF THE WAVE IN GHZ 

C SIGM IS THE STANDARD DEVIATION OF THE CANTING 

C ANGLE OF THE RAIN DROPS 
C THE IS THE AVERAGE CANTING ANGLE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcX 

COMPLEX F(2.2),EIGV(2),EICVC<2,2),WK(2),DIAG(2,2),UNIT(2,2) 
C,WA(2),EI(2, 1 ),EFIN(2,1 ) , EGVC( 2, 2 ) , MATR1 ( 2, 2 ) ,MATR2( 2, 2 ) 

C, MATR3 (2,2) 

COMPLEX FH. FV. FSII, JAI , SCI HE, Dl FF, ER1 , I SOL, ATT 
C, ERO, ELO, FV1 , FV2, FH1 , FII2, D I FF1, Dl FF2, FSH1 , FSH2 
C, M I NUSJ , EXC, EYC, EXX, EYX 



cn cun c~3 


r 
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r~\ cn nzi cm m cm m a 


REAL NO, KO, LOKANG, L, LAM, ISOLDS, I SD( 52 ) , RATE( 52 ) 
C, AMODE( 6 ) , ATTD( 52 ) 

DATA AMODE/0. 0, 0 . 25, 1 . 0. 2. 0, 2. 5, 3 . 5/ 

COMMON/VAR/ I SD, RATE, AT FD 
COMMON/BLOC 1/ MINUSJ,CONV 
COMMON/ BL0C2/ VWACM 
M I NUSJ- ( 0 . 0 , - 1 . 0 ) 

P 1 = 3 . 1415927 
CONV=P I / 1 00 . 0 
J A I = ( 0, 0, 1.0) 

P1=0.9 
P2=0. 6 
C 

C READ IN THE VARIABLES OF T *? '*.« CELL 
C 

15 RFADj S,999» I OKANC, THE1 , C !unl , NO, FREQ 
I T( LOKANC. EQ.O. ) GO TO 86 
READ( 5,977) FPSW. TAUW, EPSC, TAUC, EPSX, TAUX 
977 F0RMA1 ( 6F10 . 0 ) 

V/RITE<6.93) 

WR I TE( 6, 92 ) TIIE1, SICM1.L 
C 

C READ THE INCIDENT VECTOR FIELD 









c 


C READ( 5,91) (EI(J,1),J=1,2), ISTAT 
C 91 F0RMAT(i4(F»l.1,1X).1X, 11) 

LOKANG= LOKmNG»P I / 1 80 . 



TI1E=THE1*P 1/100. 

S I CM- S I CM 1 * P I / 1 80 . 

LAM=0. 3/FREQ 
KO=2.»>PI/LAM 

CALL C0MPNT( EPSW, TAUW, E I (1 , 1 > , E I ( 2, 1 ) ) 

CALL COMPNT ( EPSC, TAUC, EXC, EYC ) 

CALL COMPN1 ( EPSX, TAUX, EXX, EYX ) 

EYC=CONJG(EYC) 

EYX-CONJG( EYX) 

VWACM=CABS( EXC*E l{1 # 1 ) +EYC*E 1(2,1)) 

CALL ERROR ( S I CM, ER1 , 2. ) 

SCTHE=CEXP( J A I *2 . *THE ) *ER1 »EXP( -2.*(SI CM**2 ) )/2. 

C 

C CAICULATE THE AVERAGE OF COS(2.G*TIIETA) , SIN(2.*THETA) 
C 




I .—I 


CZZ3 lT 





► 


* 



CTHE=REAL( SC I HE ) 

">THE=A I MAG( SCTHE ) 

STHE=0. 0 

CTHE=EXP( -2. +( SI6M* # 2 ) ) 

WRITE(6,777 ) CTHE 
777 FORMAT ( IX, FI 5. 5) 

I K=9 

DO 50 I JK=1 , 50 
R0=10. +145 . *( IJK-1 )/50. 

I K= I K+1 
DO 66 1=1,2 
DO 66 J= 1 , 2 
66 F ( l,J)=(0.0,0.0) 

DO 65 J=1 , 2 
AL=0.8*L 

R=RO*( (R0/10. )**(-0.66) ) 

I F( J . EQ . 2 ) AL=L*0. 2 
I F( J . EQ. 2 ) R=R0 
DO 65 1=1,5 
AM0DE1=AM0DE( I ) +0.001 
AM0DE2= AMODE ( 1 + 1 ) 

CAI L COFF( TREQ, AMODE 1 , FV1 , Fill , FSH1 , D I FF1 , LOKANG, R ) 





LO 

Ln 

I 


CALL COEF( FREQ, AM0DE2, FV2, FH2, FSH2, Dl FF2, LOKANG,R) 

FV=FV2-FV1 

FH=FM2-FH1 

FSH=FSH2-FSH! 

Dl FF=l)l FF2-DIFF1 

F( 1 , 1 )=F( 1 . 1 )+( P1*FSH+P2 # ( ( FH+FV)/2. -Dl FF 
C*CTHE/2. ) )*AL*NO 

F( 2,2 )=F( 2, 2 )+( P1*FSH+P2*( ( FV+FH)/2.+DI FF 
C*CTHE/2. ) )»AL*NO 

F( 1 ,2 )=F( 1 ,2)+( P2*DIFF )*AL*N0*STHE/2. 

F(2, 1 )=F» 1,2) 

65 CONTINUE 

F(1»1)=2.*PI*F(1,1)/K0 
F(2,2)=2.«PI*F(2,2)/KO 
F( 2, 1 )=2. *P I *F( 2, 1 )/KO 
F(1,2)=F(2,t) 

C 

C FIND THE EIGENVALUES AND EIGENVECTORS OF THEPROPAGAT I ON 


C TENSOR F 
C 


CZ3 CZH 


, , 




- v 







ORIGINAL PAGE 13 

OF POOR QUALITY 


WM ; i i mn t. — i 




CALL E IGCC( F,2,2, 1, E I GV, E I GVC, 2, WK. IER1 ) 

D I AG( 1,1 )=CEXP( -JAI *E I GV( 1 ) ) 

C DIAG(1,1)=1.-JAI*EIGV(1) 

C DIAG(2,2)=1 . -JAI # EIGV(2) 

01 AG( 2,2 )=CEXP( -JAI #E IGV( 2 ) ) 

DO 10 1=1,2 
DO 10 J= 1 , 2 
I F( I .EQ.J) GO TO 10 
0 1 AC ( l,J)=(0.0,0.0) 

10 CONTINUE 
00 20 1 = 1,2 
DO 20 J=1,2 
UN I T( l,J)=(0.0,0.0) 

IF(I.EQ.J) UNIT< l,J)=(1.0,0.0) 

20 CONTINUE 
DO 30 1=1,2 
DO 30 *1=1 ,2 

30 EGVC( I , J )=EIGVC( l,J) 

C 

C FIND THE INVERSE OF THE MATRIX OF EIGENVECTORS 
C AND THEN CALCULATE EXP(F*L) * INCIDENT VECTOR FIELD 

C 




CALL LEQT1C( E I GVC, 2, 2, UN I T, 2, 2, 0, WA, I ER) 

CALL MULT( EGVC, 2, 2, D I AG, 2, MATR1 ) 

CALL MULT( MATR1 , 2,2, UN I T, 2, MATR3 * 

CALL MULT ( MATR3 ,2, 2, El, 1, EFIM) 

I F( IK.LT. 10) CO TO 150 
WRITE! 6, 97) RO, FREQ 

97 FORMAT! IX,/, ' FOR RAIN RATE* , F5. 1, 'MM /HR AND FREQUENCY' , F5. 1 , 
C' GHZ THE UNDEPOLARIZED POLARIZATION DIRECTIONS ARE:') 

DO 55 J=1,2 

WR I TE( 6, 96 ) (UNIT! I, J), 1=1,2) 

96 FORMAT! IX, ' ( ' , F5.2, '+JA *,F5.2,') XO + (',F5.2, 

C'+JAI ' ,F5.2, ' ) YO' ) 

55 CONTINUE 
I K=0 

150 CONTINUE 

CALL OUT ANT ( EFIN(1,1),EFIN(2, 1 ) , EXC, EYC, EXX, EYX, I SOLDB, ATTDB 
C, PIIAS) 

I SD( I JK)= ISOLDB 
ATTO( I JK )=ATi DD 
RATE! I JK )=R0 
50 CONTINUE 


CALL MAP 


WR I T E ( 6 , 8 ) { I SD( J ) , ATTD( J ) , RATE( J ) , J=1 , 20 ) 

8 FORMAT ( 5X, 3 F 1 0 . 3 ) 

CO TO 15 

999 FORMAT ( 3 F5 . 2, F8.2, F12. 3, F6. 3 ) 

92 FORMAT ( 5X, ' THE AVERAGE CANTING ANGLE: ' , F5 . 2, 5X. 

C' THE STANDARD DEV. OF THE CANTING ANGLE: ' , F5 . 2, /, 

C5X. >HE EFFECTIVE LENGTH 1 , F8.2) 

93 FORMAT( 1H0, 15X, ' THE VARA I BLESS OF THE RAIN CLOUD ARE'//) 
88 STOP 


C 

C 

C 

C 

C 


END 


SUBROUTINE MULT( A, N, M, B, Ml , C ) 

THIS SUBROUTINE MULTIPLIES TWO COMPLEX MATRICES 


COMPLEX A(N,M),B(M,M1 ) , C( N , Ml ) , SUM 


SUM=(0. 0,0.0) 

00 20 IJ=1,M 

20 SUM=SUM+A( I, IJ)*B( IJ,J) 

10 C(I,J)=SUM 
RETURN 
END 

C 

C 

SUBROUTINE ERROR( SIGM, ER, ALFA) 

C 

C MIIS SUBROUTINE CALCULATES THE ERROR FUNCTION OF 
C A COMPLEX VARIABLE AND IT USES THESE RESULTS TO 
C CALCULATE THE COMPLEX COEFFICIENTS OF THE AVERAGE 
C COS( ALFA* PHI) AND SIN( ALFA* Pill) WHERE PHI IS 
C A RANDOM VARIABLE WITH NORMAL DISTRIBUTION (O.SIGM) 
C 

COMPLEX ER,JAI,Z1,Z2,Z,ZS,W1,W2,X3 
JAI=(0. 0,1.0) 

Z1=3 . 3553-JA I *ALFA*S I GM/SQRT( 2. ) 

Z2=-3 . 3553-JAI *Al FA*SIGM/SQRT(2. ) 

Z=JAI*Z1 
Z S=Z*Z • 
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W1=-Z1*((0.4613135/(ZS-0.1901635)) + (0.09999216/(ZS- 
C1 . 78*14927 ) )+( 0 . 002883B94/! ZS-5 . 5253437 ) ) ) 

W2=2 . *CEXP( J A I *5 . *ALFA#S I GM ) 

X3=JA I *A I MAG( Z2*Z2 ) 

ER=W2*CEXP{ -X3 )-C0NJG(W1 )*CEXP( -Z2*Z2 )-W1*CEXP( -Z1*Z1 ) 
RETURN 
END 
C 
C 

, SUBROUTINE MAP 

£ C 

? 

C THIS SUBROUTINE PI OTS ATTENUATION AND 
C ISOLATION VERSUS RAIN RATE 

C 

REAL I SD( 52 ) , RATE( 52 ) , ATTD( 52) 

COMMON/VAR/ I SD, RATE, ATTD 
CALL PLOTS( 0,0, 50 ) 

CALL SCALE ( RATE, 8.0,50, 1 ) 

»T? - t { ;•*_ f ... t . * . -*w ' ?> 

CALL SCALE! ISO, 6. 0,50, 1 ) 

CALL SCALE! A1TD, 6. 0,50, 1 ) 

CALL PLOT! 2. 0,2.0, -3) 

CALL AXIS( 0.0,0.0, 9IIRAI N RATE,-9,8.0,0.0,RATE(51 ),RATE(52) ) 


S o 

O 


f 2 


• ■ ■ y* -A •»# *< 




- 141 - 



CALL L I NE( RATE, I SD, 50, 1,0,0) 


CALL PLOT (0.0, 0.(1, +999) 


CALL PLOTS(0,0,50) 


CALL PL0I(2. 0,2. 0,-3) 


CALL AX I S( 0.0, 0. 0, 9HRAI N RATE, -9, 8 .0,0. 0, RATE( 51 ) , RATE( 52 ) ) 


CALL AX IS(0. 0,0.0, 1 1 IIATTENUAT I ON, 11, 6.0, 90. 0, ATTD( 51 ) , ATTD( 52 ) ) 


CALL LINE(RATE,ATTD, 50, 1,0,0) 


CALL PLOT (0.0, 0.0, +999) 


RETURN 


SUUROUT I NE COEF (FREQ, AMODE, FV, FH, FSPH, D I FF, LOKANG, R ) 


THIS SUBROUTINE RETURNS TO THE CALLING PROGRAM THE SCATTERING 


COEFFICIENTS FOR SPHERICAL AND OBLATE RAIN DROPS. 


THE COEFFICIENTS ARE A FUNCTION OF FREQUENCY AND DROP SIZE 


AND ELEVATION ANGLE. THE COEFFICIENTS USED ARE THOSE OF UZUNOGLU, 


EVANS AND HOLT. THIS ISA MODIFIED VERSION OF THE 


*.V /sir-.’ 


U 

O 5 

73 xr 


iO 73 

P fR 

Zt ip 




rm 



C SUBROUTINE DEVELOPED BY PRESINGER AND STUT2MAN. 

C 

COMPLEX CMPLX 
REAL IOKANG 
C 

COMPLEX FV.FH.FSPH.DIFF 

DOUBLE PRECISION U22, U33, U44, U55, AK1 , AK2, AK3 , AK4, AK5 
AK=-8.2*( R##( -0.21 ) ) 

AK1=EXP( AK*AMODE )/AK 
C 

I F( AMODE.LT. 0.25) AM0DE=0. 25 
I F ( AM0DC.LE.0.25) GO TO 100 

I F( ( I NT C FRFQ ) . EQ. 1 1 ) . AND . ( AMODE . GT. 3 . 5 ) ) AM0DE=3 . 5 
I F ( ( I NT ( FREQ ) . EQ . 1 4 ) . AND . ( AMODE . GT . 3 . 5 ) ) AM0DE=3 . 5 
I F ( ( I NT ( FREQ). EQ. 20). AND. (AMODE. GT. 3.0)) AMODE=3.0 
I F ( ( I NT( FREQ ) . EQ . 30 ) . AND . { AMODE . GT . 3 . 0 ) ) AM0DE=3 . 0 
IF(( INT(FREQ).EQ. 20). AND. (AMODE. GT. 3.0)) GO TO 100 
IF( ( INT( FREQ). EQ. 30). AND. (AMODE. GT. 3.0)) GO TO 100 
C 

U1 1= AMODE 
U22=AM0DE**2 
U33=AMODE**3 





'• ... » nllhaiM ji * 4 




U4M=AM0DE**4 

U55=AMODE**5 

AK2=AK**2 

AK3=AK # * 3 

AK4=AK*#4 

AK5=AK**5 

U1=( U11-1 ,/AK) 

U2=( U22-2 .*U11 /AK+2 . 0/AK2 ) 

U3= ( U33- 3 . “U22/AK+6 . *U 1 1 /AK2 
C-6.0/AK3) 

Ul|=( U44-4. *U33/AK+12. *U22/AK2- 
C2«l.*Un/AK3+2l|./AKI|) 

U5=( U55-5. *U4»l/AK+20. *U33/AK2 
C-60. *U22/AK3+120. *111 1/AK4- 120 . /AK5 ) 



GO TO 101 
ICO U1=AM()DE 


U2=(AMODE**2) 
U3=( AM0DE**3 ) 
UI|=(AMOOE*#4) 





era m a m m rm 



i 

h-* 




i 




U5=(AMODE*#5) 

101 CONTINUE 
C 
C 
C 

C 11.0 GHZ COEFFICIENTS 

C 

I F ( I NT ( FREQ) . NE. 1 1 ) GO TO 1 

C 

C SPHERICAL DROP COEFFICIENTS 

C 

I F( AMODE.GI. 1 .00) GO TO 10 

C 

E0R=-0. 002054 8+0. 01 638947*U 1-0. 04 1 7568 # U2+0. 0883221 3*U 3 
F0f=-0. 0025154+0. 01928553*U1-0.0456816*U2+0.03655147*U3 
GO TO 11 
10 CONTINUE 
C 

E0R=-1. 281 55706+2. 83287718*U1 -2. 07399678»U2+0. 60 1.90887*U3 
1-0. 00984 194*um-0.01096165#U5 
£01=2.60278025-7. 52434662*U1+8. 14691632*U2-4 . 121 33971 #U3 
1+0. 9908946 /*U4-0. 0873 1788#U5 
C 




GO TO 12 


11 CONTINUE 
C 

C OBLATE DROP COEFFICIENTS 
C 

EV90R=-0. 001 322+0. 01 036867*U1 -0. 025372*U2+0. 0707253 3 *U3 

EV90 I =-0. 0024306+0. 0186 1967*U1 -0. 0440232*U2+0. 0350741 3 *U3 

LII90R=-0.0023684+0. 01878307*01 -0.0473704*U2+0. 09255573*03 

EH90l=-0. 0030366+0. 02320527*U1-0.0546296*U2+0.04299093*U3 

0 1 1 FR=-0 . 0010464+0. 0084 144*vH -0 . 0219984*U2+0. 02 1 8304 *U3 

£ D I FF I =-0 . 000606+0 . 0045856*U 1 -0 .0106064*02+0 . 0079168*03 

Ln 

1 C 

GO TO 1000 

12 CONTINUE 
C 

EV90R=-0. 3653892+0. 32540943*01+0. 45280976*U2-0. 5^907398*03 
1+0. 2030921 5*U4-0.02432204*U5 

EV90 1=2. 20618555-6. 19251766*U1+6.43873903*U2-3.090414*U3 
1+0. 7060021 1*U4-0. 0607 173*U5 

EH90R=2. 09022 F06-7. 073451 19*U1 +9. 03639726*U2-5. 33 15621 5*U3 
1+1 .49Yl6423*U4-0. 15747063*05 

FH90 1=3. 28654422-9. 682 1 7951*01+10. 71102231*U2-5. 5491 184*U3 
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c 

c 

c 

c 

c 


GO TO 1000 
i CONTINUi 

1H GHZ GOEFF ICI ENTS 

















era □ a czn n i — i 






■sfe I 





EV90R=-0 . 0002*18+0 . 00335967*01 -0. 01 3928*02+0 . 10002 1 33*U3 
EV90 I =-0. 008936+0. 06*155067*1)1-0. 152256*02+0. 1 1893333*03 
EH90R=-0. 000656+0. 007656*01 -0, 029656*02+0. 122656*03 
EII901=-0. 010366+0. 079159*01-0. 18596*02+0, 193872*03 
01 1 f R--0. 000908+0. 00930 133*01-0.01 5728*02+0.02263967*03 
D I TF I =-0.001 93+0. 01 959533*01 -0.033709*02+0. 02993867*03 
C 

GO TO 1000 
22 CONTI NOE 
C 

EV90R=-9. 98993663+12. 95981829*01 -12. 90182271 *02+6. 3290 191 3*03 

1-1.92805671*09+0. 1 1875301*05 

EV90l=-0. 1923376+0. 3925132*01-0. 196292*02+0.07968968*03 
C 

I F ( AMODE . GT. 2. 5 ) GO TO 29 
C 

E090 1=1.1 0278958-3 . 79696237*01+9 . 98092001*02-2.2296931 7*03 
1+0.92010037*09 

01 FFR=-0. 15268909+0. 39830895*01-0. 3527577*02+0. 1 1767019*03 
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D i FF I =- 1 . 2996004 1 +2 . 93303524*U1 -2 . 1 5255009*U2+0. 5206384*U3 


I F ( AMODE . GT . 2 . 00 ) CO TO 25 
C 

EH90P.= -1 .02599992+2. 307333 17*U1-1 . 60799988*U2+0. 4266664»U3 
C 

GO TO 1000 

24 CONTINUE 
C 

EII90 1=2.51 000033-2 . 03000022*U1 +0 . 62000004*U2 
Dl TFR=-15. 85152202+18. 27189723*U1-6.68140666 # U2+0. 7941 1 153*U3 
01 FFI=22. 32321586-44. 48681458*U1+30. 17926l25*U2-8. 46925565*U3 
1+0. 85300427*U4 
C 

CO TO 1000 

25 CONTINUE 

C 

EH90R=- 19 . 65+20. 6567* U1 -6 . 66*U2+0. 6933*U3 
C 

GO TO 1000 
2 CONTINUE 
C 

C 20 GIIZ COEFFICIENTS 
C 


O O 
-n 33 

-O § 
O 2 

O V-- 

73 r* 

o -0 
c 

> © 
r 


^ ■■ 


dll 


L=Z3 




i 




— , — — - — ■ - - ■ ' r 
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I F( I NT( FREQ) . NE. 20 ) GO TO 3 
C 

C SPHERICAL DROP COEFFICIENTS 
C 

I F( AMODE.GT. 1.0) GO TO 30 
C 

E0R=0 . 020296-0 . 1 45276*U 1 +0 . 297656#U2+0 . 008224*U3 
C 

EO I =-0.0 159 88+0. 12709133*U1-0. 334568*U2+0. 30563467*03 
C 

GO TO 31 

30 CONTINUE 
C 

EOR= 3.355671 52-8 .6181 8659*01+7 . 8 1 57909*U2-2 . 682 1 79 1 7*U3 
1+0. 3108073*U4 

EO 1=1 .85636463-3. 70750898*Ul+2. 19599302*U2-0. 25148145*U3 
C 

CO TO 32 

31 CONTINUE 

C 

C OBLATE COEFFICIENTS 


C 
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EV90R=0. 02454800-0. 1 79121 33*01+0.384448*02-0.07447467*03 
EV90l=-0. 013492+0. 11190533*01-0.299032*02+0.27853867*03 
EII90R=0. 02282-0. 16596667*01+0. 34976*02-0.03061333*03 
EH90l=-0. 015178+0. 12623933*01-0. 33848*02+0. 31533867*03 
0 1 FFR=0 . 000072+0 . 001 75467*0 1 -0. 01 5488*02+0 . 03426133*03 
0 IFF l=-0. 001686+0. 014334*U1 -0. 039448*02+0. 0368*U3 

GO TO 1000 
32 CONTI NOE 


EV90R=1. 56037246-3. 97019757*01+3. 66524972*02-1. 24458239*03 
1+0. 14526662*04 

EV90 1=4. 77376463-13.63610787*01+14.83397997*02 
1-7.66136891*03+1.96881573*04-0. 19749381*05 
EH90R=4. 63026679-1 1 .97374109*01+1 1.0254653*02 
1 -4 . 00506 13*03+0.50244373*04 

EII90I=- 1.5390057+6. 3417866*01-9. 35328879*02+6. 20987895*03 
1-1 . 74544867*04*0. 17820925*05 
D I FFR=3 . 06989428-8 . 00354339*01+7 . 36021547*02 




rrm mm cm rrm cm im im mm 


mm m cm 


m BSB® 


1-2.76047886»U3+0. 3571771*U4 

Dl FF l=-6. 3 1274803+19. 97 7833 34 # U1 -24 . 187205 3*U2+1 3 . 8712166*1)3 
1-3.71425705*U4+0.3757024*U5 
C 

GO TO 1000 
3 CONTINUE 
C 

I F ( INT( FREQ). NE. 30) GO TO 2000 
C 

C 30 GHZ COEFFICIENTS 

l 

(— 1 r 

Ui C 

PO 

1 C SPHERICAL DROP COEFFICIENTS 

C 

I F ( AMODE . GT. 1 . 0 ) GO TO 40 
C 

E0R=0. 028004-0. 221 564+U1+0. 530744*U2+0.022816*U3 
E0l=-0. 007072+0. 09329866*U1 -0. 3901 12#U2+0. 55<>98533*U3 
C 

GO TO 41 
40 CONTINUE 
C 


E0R=-1 .95096901+3. 4107597*U1-0.82790307*U2-0. 39982485*U3 
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1+0. 1291964*U4 

E0I=5. 77798066-16. 3876166*01+16. 10102407*02-6. 06482662*113 
1+0. 82747918*114 


Ui 

LO 

I 


GO TO 42 
41 CONTINUE 

OBLATE DROP COEFFICIENTS 

EV90R=0. 03542-0. 2838733*01+0. 7038*02-0. 15834666*03 
EV90 1=0. 01 3876-0. 064756*U1-0.026904*U2+0. 290784*03 
EII90R=0 . 03 1 652-0 . 252652*111 +0 . 6 14952*U2-0 . 049952*113 
EII90l=-0. 000532+0. 04746533*01-0. 299512*U2+0.50957867*U3 
DIFFR=-0. 003768+0.03122133*01-0. 088848*112+0. 10839467*03 
01 FFI--0. 014408+0. 1 12221 33*01-0. 272608*U2+0. 2187946* ***'. 


O O 
■n X 

Tj G 

o9 

O 

to r 
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c > 
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GO TO 1000 
42 CONTINUE 


EV90R=3 . 93 161268- 10. 97208846*U1+1 1 . 43870068*U2 
1-4 . 78953064*03+0. 69447892*114 
EV90I=3. 12729055-9. 08249632*U1+9. 16246329*U2 


cm a cm m 


m m 


rm rm 


1-3.47227982*03+0.47789632*04 



i 


Ul 

I 


EH90R=-6 . 97583584 + 16 . 79179963*01-13.53607686*02 
1+4.6381 185l'*J3-0. 57535653*04 

EH90I=-1. 39373246+1. 72736066*01-0. 14076091*02+0.05474057*03 
01 FFR=-1 0.9074465 1+27. 76388807*01-24.97477753*02 
1+9. 42764918*03-1 . 26983545*04 
0 1 FF I =5 . 33984552- 1 7 . 20920477*01+20 . 781254 1 *02 
1-11. 64360435*03+3.07942984*04-0. 30370378*05 
C 

1000 CONTINOE 

€ 

ALPII=1 . 576796-LOKANG 
CSLA=0.001*C0S( ALPII J**2 
SNEA=0.001*SIN(ALPH)**2 
FVR=CSLA*E0R+SNLA*EV90R 
FVI =-CSI A*E0 1 -SNLA*EV90 1 
FHR=CSLA*E0R+SNLA*EH90R 
FH I =-CSI.A*E0 I -SNI.A*EH90 1 
DirFR--SNLA*DIFFR 


If 

|5 






j ; UbstLi&xmAiL. • 


. 






DIFFI=SNLA*DIFFI 

E0R=0.001*E0R 

E0I=-0.001*E0I 

C 

FV-CMPLX( FVR, FVI )*AK1 
FH=CMPLX( FHR, FHI )*AK1 
FSPH=CMPLX( EOR, EO I )*AK1 
Dl FF=CMPLX(OI FFR.DI FFI )*AK1 
C 

GO TO 3000 

Go. 

Ln 

‘ 2000 WRITE(6,2001 ) 

2001 F0RMAT(//,3X, 'FREQUENCY NOT ALLOWED, ONLY 11,14,20,30 GHZ ALLOWED' 

1) 

C 

STOP 

3000 CONTINUE 
C 

RETURN 

END 

C 

C 


c=3 a Em □ c=3 ::j :m c=! i □ a~ tr cn cm cm era =3 n=a t=n 
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SUBROUTINE OUTANT ( EXN, EYN, EXC, EYC, EXX, EYX, I SOLI , ATTEN1 , PHASE1 ) 

THIS SUBROUTINE TAKES THE X AND Y COMPONENTS OF THE WAVE EXITING 
THE RAIN CELL (EXN.EYN) AND USES THE COMPLEX VECTOR METHOD 
TO COMPUTE VALUES FOR ATTENUATION, ISOLATION, AND PHASE AS A 
RESULT OF THE RAIN MEDIUM AND POLARIZATION MISMATCH EFFECTS OF THE 
RECEIVE ANTENNA 

THIS DATA IS STORED IN PROGRAM MEMORY FOR LATER OUTPUT 
THIS SUBROUTINE WAS WRITTEN BY PRESINGER. 

REAL I SOLI 

COMMON/BLOC 1/MINUSJ, CONV 
COMMON/BLOC2/ VWACM 
COMPLEX MINUSJ 
COMPLEX EXC, EYC, EXX, EYX 


COMPLEX VWPAC, VWPAX 
COMPLEX EXN.EYN 
REAL VWPACM, VWPAXM 







ft* 



Ln 

I 


VWPAC=EXN*EXC + EYN*EYC 
VWPAX=EXN*EXX + EYN*EYX 


VW P ACR= R E A L ( VW P AC ) 

VWPAC I = A I MAG( VWPAC ) 

VWPACM=CABS( VWPAC ) 

VWPACR= ATAN2 ( VWPAC I , VWPACR ) 

VWPAC P=VWPACR/CONV 

VWPAXR=REAL( VWPAX) 

VWPAX I = A I MAG( VWPAX ) 

VWPAXM=CABS( VWPAX ) 

I F( VWPAXM. EQ. 0. 0 ) GO TO 10 
VWPAXR=ATAN2( VWPAX I ,VWPAXR) 

VWPAX P=VWPAXR/CONV 
I F( VWPAXM. LE.O. 000001 ) VWPAXP=0.0 


T) § 

85; 

to r 

<o 3 

c r* 

3 .Z 1 


I SOL 1 =20 . w ALOG 1 0( VWPACM/VWPAXM ) 
GO TO 11 
10 CONTINUE 
VWPAXP=0.0 


cm cm m 


t=a 


-8SI- 


rm r— j 


m 


m r~3 m n 


I S0L1=999. 99 
11 CONTINUE 
C 

PHASE1-VWPAXP-VWPACP 

C 

I F( PHASE 1 . LT . 0 . 0 ) PHASE 1=PHASE1 +360 . 0 
C 

AT TEN 1 =20 . *AL0G10( VWACM/VWPACM ) 

C 

RETURN 

END 

C 

C 

SUDROUT I NE COMPNT( EPS, TAU, EX, EY) 

C 

C THIS SUBROUTINE RETURNS THE X AND Y COMPONENTS GIVEN AN EPSILON 
C ANO TAU (IN DEGREES) DESCRIBING AN ARBITRARY POLARIZATION STATE 
C 

COMPLEX EX,EY 
COMPLEX MINUS.) 

COMMON/BLOC 1 /M I NUSJ , CONV 

I 

C 


O o 
^ XJ 

r? & 

O s 
-o r . 
'to 

£ Q 

C tr { 
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EPSR=EPS*CONV 
TAUR=T AU*CONV 

IF(ABS(EPSR).EQ.(45.*CONV)) GO TO 1 
I F( EPSR. EQ.O. ) GO TO 2 
1 F( TAOR. EQ.O. > GO TO 3 
I F ( TAUR. EQ. ( 90. *CONV) ) GO TO 4 



T 1=1 AN ( 2 . *EPSR ) 

T2=SIN(2.*TAUR) . 

DEl.TR=ATAN2(T1,T2) 

CAMR=0 . 5#ARC0S( COS( 2*EPSR )#COS( 2*TA0R ) ) 
GO TO 100 

1 OEl.TR-2. *EPSR 
GAMR=45. M CONV 
GO TO 100 

2 OELTR- 0 . 

GAMR=1AUR 
GO TO 100 


g=i =i =a 


=□ 


1 
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3 DEI..TR=S!GN( 1 . , EPSR) *90. *CONV 
GAMR=ABS( EPSR ) 

GO TO 100 

»• DCLTR=SIGN( 1 . , EPSR )*90. *CONV 
GAMR=90 . *CONV-ABS( EPSR ) 

100 CONTINUE 
C 

EX=COS( GAMR ) 

EY=S I N( CAMR ) *CEXP( -M I NUSJ*D£LTR ) 
C 

RETURN 

END 



